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Multi-Objective Genetic Algorithm to Design 
Manufacturing Process Line Including Feasible and 
Infeasible Solutions in Neighborhood 

Masahiro Arakawa and Takumi Wada 

Department of Civil Engineering and Systems Management, Nagoya Institute of Technology, Gokiso, Showa, Nagoya, 466-8555, 
Aichi, Japan 
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Abstract: This paper treats multi-objective problem for manufacturing process design. A purpose of the process design is to decide 
combinations of work elements assigned to different work centers. Multiple work elements are ordinarily assigned to each center. 
Here, infeasible solutions are easily generated by precedence relationship of work elements in process design. The number of 
infeasible solutions generated is ordinarily larger than that of feasible solutions generated in the process. Therefore, feasible and 
infeasible solutions are located in any neighborhood in solution space. It is difficult to seek high quality Pareto solutions in this 
problem by using conventional multi-objective evolutional algorithms. We consider that the problem includes difficulty to seek high 
quality solutions by the following characteristics: (1) Since infeasible solutions are resemble to good feasible solutions, many 
infeasible solutions which have good values of objective functions are easily sought in the search process, (2) Infeasible solutions are 
useful to select new variable conditions generating good feasible solutions in search process. In this study, a multi-objective genetic 
algorithm including local search is proposed using these characteristics. Maximum value of average operation times and maximum 
value of dispersion of operation time in all work centers are used as objective functions to promote productivity. The optimal 
weighted coefficient is introduced to control the ratio of feasible solutions to all solutions selected in crossover and selection process 
in the algorithm. This paper shows the effectiveness of the proposed algorithm on simple model. 

Keywords: Process design, process line, feasible and infeasible solution, multi-objective genetic algorithm, mix production, 
simulation 



1. Introduction 

This paper treats multi-objective problem for 
process design. A purpose of process design is to 
decide combinations of work elements assigned to all 
work centers. Here, the work elements denote works 
to produce a single product and multiple work 
elements are ordinarily assigned to every center. There 
is possibility that infeasible solutions are generated by 
precedence relationship of work elements in process 
design. In the process of process design, the number 
of infeasible solutions sought is ordinarily larger than 
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that of feasible solutions. However, there are the cases 
that infeasible solutions are superior to feasible 
solutions in terms of objective functions. 

This problem is treated as multi-objective problem 
in order to design process line for high productivity. 
Especially mix production appeared in a large number 
of process lines is focused. The following functions 
are introduced as multi-objective functions to evaluate 
the productivity of process line for mix production: 
the maximum value of average operation times in all 
work centers and the maximum value of dispersion of 
operation times in all work centers. 

However, this problem is including feasible and 
infeasible solutions in any neighborhood in solution 
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space. Therefore, it is difficult to seek high quality 
feasible solutions by using evolutional algorithms. 

In past studies related to process design, many 
researchers proposed heuristic methods and rule -based 
procedures [1-6]. Here, one of the most conventional 
and popular procedures for process design is the 
ranked positional weight technique [6]. This is a 
heuristic procedure that has been used to develop 
many designs. However, process design is focused on 
a single-objective problem for line balance in these 
studies. In addition, infeasible solutions are neglected 
or infeasible solutions are modified to feasible 
solutions in the search process. In the case of 
conventional genetic algorithm, lethal chromosomes 
are modified to normal chromosomes to prevent from 
generating infeasible solutions. 

We consider that a problem of process design 
includes the following characteristic causing difficulty 
to seek high quality solutions: 

(1) Many infeasible solutions which have good 
values of objective functions are easily sought in the 
search process, since infeasible solutions are resemble 
to good feasible solutions, 

(2) Infeasible solutions are useful to select new 
variable conditions generating good feasible solutions 
in search process. 

When new conditions are generated from existing 
solutions in evolutional algorithm under the condition 
that feasible and infeasible solutions have same 
possibility to be selected, Characteristic (1) causes 
generation of a large number of infeasible solutions in 
the search process. 

On the other hand, a genetic operation using 
Characteristic (2) is expected to promote seeking high 
quality feasible solutions since various patterns of 
chromosomes are generated. Therefore, 
multi-objective genetic algorithm (MOGA) including 
a search procedure using Characteristic (2) is 
proposed in this paper. Both chromosomes which 
generate feasible and infeasible solutions are used to 
select as parent chromosomes in crossover process 



and as chromosomes preserved for next generation in 
the proposed algorithm. Probability to select 
chromosomes generating infeasible solutions is 
controlled in crossover and selection for next 
generation in the proposed algorithm. Then, high 
quality feasible solutions are expected to be sought. 

In this paper, the characteristics of the proposed 
algorithm are explained and the performance of the 
proposed algorithm is evaluated on simple model for 
process design for mix production. 

2. Process Design 

2.1 Characteristics of Process Design and Objective 
Functions 

Fig.2 shows schematic diagram of assignment 
process for process design. Process design denotes 
that all work elements required to manufacture 
products are assigned to work centers. A precedence 
relationship is constructed among different work 
elements and this relationship is affected as 
constraints to assign the work elements to the work 
centers. Therefore, when a procedure relationship is 
not satisfied between different work elements assigned 
to works, infeasible solutions are generated. 
Ordinarily, since multiple work elements are assigned 
to each single work center, infeasible solutions which 
take good resultant values of objective functions are 
easily generated. 

O O Infeasible 
solution 

Feasible 
solution 



£2 




fl 

Fig. 1 Schematic diagram of search process for process 
design. 
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Operation 
Time 



Practical Cycle Time 




Total operation time to produce a product 

Fig. 2 Schematic diagram of assignment process for 
process design. 

Minimum of work centers is calculated from 
production plan. In conventional process design, 
target cycle time is calculated to satisfy production 
plan. Here, target cycle time denotes time interval to 
manufacture required products from a production line 
in predetermined period. The cycle time is defined as 
the maximum operation time of work elements 
assigned among work centers. 

The minimum number of work centers 'Nmin' and 
target cycle time 'Tc' are calculated from Eqs. (1) and 
(2) from production plan. 



N„ 



\TJT C 



T c = S/Q 



(1) 

(2) 



Here, T all is total operation time to complete a 
product. S is planning period where products are able 
to be produced. Q is total amount of products required 

to produce. [*~] denotes floor function. 

A purpose of conventional process design is to 
assign work elements to work centers under the 
condition that maximum operation time among work 
centers is equal to or smaller than the target cycle time. 
Process line calculated under the condition easily 
causes large difference of operation times among 
work centers. Since maximum operation time among 
work centers is practical cycle time, minimizing the 
practical cycle time is required to promote 
productivity. Large variation of operation time causes 
idle time at continuous work centers. Therefore, 
minimizing variation of operation time at all work 
centers is effective not to decrease productivity. 

From this discussion, maximum of operation time 



in work centers 'Max(/i n y and dispersion of operation 
time in work centers 'Max(s w )' are introduced as 
objective functions to evaluate the quality of process 
design in this study. These functions need to minimize 
simultaneously to increase productivity. Eqs. (3) and 
(4) denote the maximum value of average operation 
time and the maximum value of dispersion of 
operation time in all work center. These are objective 
functions. 'tri denotes work center. 

Max (jj. m ) 

Max (s m ) 



Mm 



(3) 
(4) 

(5) 



Vm 



(6) 



Here, /u t and s t denote average operation time and 
dispersion of operation time of work element 7' 
respectively. d m is 0-1 variable which takes 1 when 
work element T is assigned to work center 'm.' 
Otherwise, the variable takes 0. 

In Fig. 2, arrows show direction of operation time at 
work centers to improve productivity in the process of 
assignment of work elements. When maximum of 
average operation time in all work centers is decreased 
by change of work elements among work centers, 
practical cycle time is reduced and productivity is 
increased. In addition, when maximum of dispersion 
of operation time in all work centers is decreased, 
standard deviation of operation time is decreased and 
idle time caused by variation of operation time among 
different work centers is reduced. 

The mathematical model of this problem is written 
as following: 

MinC mm (7) 
MinD^ (8) 
Subject to 

C max >is ikPi Vk (9) 



i=i 



D > 

max 



ik S i 



Vk (10) 



1=1 
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k=l 



k + M(5 ik -\)< Ul 
Ui <K(l-S ik ) + k 

U;<U; 



Mi 

Vi.VJfc 
Vi, \/k 



(11) 

(12) 
(13) 
(14) 



Here, Eqs. (7) and (8) denote objective functions. 
C max and D max denote maximum value of average 
operation time and maximum value of dispersion of 
operation time among work centers respectively, k and 
i denote work center's number and work element's 
number respectively. In this model, both k and i take 
integer numbers from one. K and / denote the 
maximum number of work centers and the maximum 
number of work elements. In addition, w, denotes a 
variable which presents work center's number to 
which work element i is assigned. M is a large number 
and it takes 2K in this problem. When this 
mathematical model is resolved, the optimal solutions 
are obtained and the solutions are used to compare 
with the proposed algorithm. 

2.2 Bi-Objective Genetic Algorithm Including 
Evaluation of Infeasible Solutions 

New type of multi-objective genetic algorithm is 
proposed to resolve bi-objective problem for process 
design in this study. In a problem of process design, it 
is known that many infeasible solutions are ordinarily 
located nearby feasible solutions in solution space 
from numerical experiment in past. Therefore, it is not 
easy to seek approximate optimal Pareto solutions by 
using conventional multi-objective genetic algorithm. 
Many infeasible solutions are ordinarily superior to 
feasible solutions in terms of resultant values of 
objective functions. The proposed algorithm uses 
infeasible solutions in order to effectively seek 
approximate optimal Pareto solutions. Both feasible 
and infeasible solutions are adopted as chromosomes 
to select for crossover process and to preserve 
chromosomes for next generation. Values of criteria of 
infeasible solutions are modified by Eq. (15) in order 



to control probability of selection of infeasible 
solutions. 

Fitness(i) = KX*value(i) (15) 
Here, Fitness(i) and value(i) denote fitness of 
chromosome 7' used for selection and resultant value 
of criterion of chromosome 7' respectively. KX 
denotes weighted coefficient to modify value of 
criterion of infeasible solutions. KX takes 1 or a 
number more than 1 . A simple MOGA starts after KX 
is determined in the proposed algorithm. By 
introducing the coefficient, we consider that an excess 
of infeasible solutions is suppressed and solutions 
which have smaller values of objective functions are 
effectively sought. The process of the proposed 
MOGA is described in the following. 

(1) Weighted coefficients located at left side and 
right side, 'KLO' and 'KRO', are determined to be 
initial values. Here, KLO is smaller than KRO. 
Parameter to control weighted coefficient number 'a' 
is determined to be constant number and the 
maximum number of iteration to change the weighted 
coefficient 'MAXITER' is determined to be a constant 
number. 

(2) The number of iteration is determined to be 0. 
The coefficient located at left side 'KL' and the 
coefficient located at right side 'KR' are calculated 
from Eqs. (16) and (17). 

KL = a (KR0-KL0)+KL0 (16) 
KR = KRO -a (KR0-KL0) (17) 

(3) The ranked positional weight technique is 
executed to generated an initial solution. 

(4) Initial population is generated by using the 
solution generate in Step (3). 

(5) The coefficient located at left side 'KL' is 
assigned to KX. Then, a single MOGA is executed 
using the coefficient determined to be KX. 

(6) The coefficient located at right side 'KR' is 
assigned to KX. Then, a single MOGA is executed 
using the coefficient determined to be KX. Move to 
Step (8). 

(7) A single MOGA is executed using the 
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coefficient determined to be KX. 

(8) Pareto solutions obtained using the coefficient 
located at left side 'KL' is compared with that 
obtained using the coefficient located at right side 
'KR.' 

(9) Move to Step (10) when feasible solutions 
obtained by using the coefficient located at left side 
'KL' are superior to feasible solutions obtained by 
using the coefficient located at right side 'KR.' Move 
to Step (11) when feasible solutions obtained by using 
the coefficient located at right side 'KR' are superior 
to solutions at the other side. Otherwise, move to Step 
(12). 

(10) Each weighted coefficient is changed by using 
the following equations: 

KRO = KR, 

KR = KL, (18) 
KL = KL0+ec(KR0 — KLO). 
Then, the coefficient located at left side 'KL' is 
assigned to KX. Move to Step (12). 

(11) Each weighted coefficient is changed by using 
the following equations: 

KL0=KL, 

KL=KR, (19) 
KR=KR0 - ec(KR0 - KLO). 
Then, the coefficient located at right side 'KR' is 
assigned to KX. 

(12) The number of iteration is added one. 

(13) Complete the algorithm when the number of 
iterations is equal to MAXITER. Otherwise, go back 
to Step (7). 

In this study, a takes 0.3820. Here, this search 
process to control one dimensional variable using this 
constant number is similar to 'Golden search method.' 
Pareto and feasible solutions obtained by using 
different coefficients are compared with each other in 
Steps (8) and (9). When the number of the Pareto and 
feasible solutions obtained by using the weighted 
coefficient KX is larger than that of the Pareto and 
feasible solutions obtained by using the other 
coefficient, the coefficient KX regards as superior 



condition to the other coefficient. 

The process of a single MOGA in Steps (5), (6), 
and (7) is described in the following. 

(51) Population 'POP' and maximum generation 
'MaxGen' are predetermined. 

(52) Current generation takes 0. 

(53) The weighted coefficient 'KX' is 
predetermined. 

(54) Crossover and mutation: 

(S4-1) Select parent chromosomes by roulette 
selection. 

(S4-2) Generate new chromosomes by crossover. 

(S4-3) Investigate whether new chromosomes are 
identical to chromosomes included in the population 
or not. If a new chromosome is NOT identical to any 
chromosome, go to Step (S5) after the chromosome is 
put into the population. When there is a chromosome 
which is identical to a new chromosome, go to Step 
(S4-4) in terms of the new chromosome. 

(S4-4) Execute mutation process for the new 
chromosome and go back to Step (S4-3). 

(55) Encode the new chromosomes and calculate 
values of objective functions of solutions obtained 
from the chromosomes. 

(56) Selection of chromosomes for new generation: 
(S6-1) Select Pareto solutions; 

(S6-2) Parallel selection; and 
(S6-3) Roulette selection. 

(57) Judge completion of algorithm. 

Complete the algorithm when the current generation 
is equal to the maximum number of generations 
'MAXGen.' Otherwise, go back to Step (S4) after one 
is added to the current generation. 

Population is continuously used under the condition 
of different weighted coefficients. Therefore, solutions 
are continuously improved in population while 
weighted coefficient is changed in the proposed 
algorithm. 

In roulette selection in Step (S6-3), values of 
objective functions of a solution obtained from a 
chromosome are modified to one dimensional 
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criterion by using Eqs.(20) and (21). Furthermore, 
inverse of the one dimensional criterion is used as 
probability to select chromosomes in roulette selection 
because the objective functions are minimized in this 
problem. 

/ 



fi 



fi+fx 

f -f 
f - f 



(20) 



(21) 



Fig. 3 shows a sample of representation of a 
chromosome. The chromosome is constructed of an 
array of work element numbers and an array of work 
center numbers. Numbers located in the same element 
number in both arrays denote that work element 
represented in the element is assigned to work center 
represented in the other element. In an array of work 
center number, numbers are arranged in the ascending 
order. In order to operate elements in an array of work 
center's numbers for crossover and mutation, new type 
of array is generated by the following modification: 1 
is assigned to the last element in each partial array of 
same work center number and 0 is assigned to other 
elements in the new array. Fig. 4 denotes an example 
of the new array. The number of '1' included in the 
array is same as the number of work centers. Therefore, 
work centers to which work elements are assigned can 
be controlled by operating locations of '1' in the array. 

In Step (S4-1), several types of criterion are used in 
roulette selection: inverse of every objective function 
and inverse of the one dimensional criterion defined 
by Eqs.(20) and (21). New chromosomes equally are 
generated in crossover which uses different criterion 
in roulette selection. When bi-objective functions in 
the problem are/! and/ 2 , > W% an d 1// calculated 
by Eqs. (20) and (21) are used as criterion in the 
roulette selection. 

In the crossover in Step (S4-2), two-point crossover 
is adopted to operate elements in an array of work 
element numbers and an array of work center numbers 



independently. Two elements are arbitrarily selected 
in the array of two chromosomes and parts of the array 
between the selected elements are exchanged with 
each other. Lethal genes are prevented by the 
following process: (a) in an array of work element 
numbers, when numbers in the original partial array 
are included in the exchanged partial array, elements 
of the numbers in the exchanged partial array are used. 
The other numbers in the original partial array are 
rearranged in the exchanged partial array according to 
the order of the numbers in the original partial array, 
(b) in an array of work center numbers, '1' is deleted at 
any elements in the exchanged partial array when the 
number of 1' in the exchanged partial array is larger 
than that in the original partial array. On the other 
hand, T is added at any element in the exchanged 
partial array when the number of T in the exchanged 
partial array is smaller than that in the original partial 
array. 

As for mutation process in Step (S4-4), swap is 
performed between any two elements by iteration of 
random number. Here, the random number is determined 
between 1 and 10 by uniform random number. 

Solutions which have different chromosome and 
same values of objective functions are stored in 
population. Population is continuously used and 
updated in all process of the proposed algorithm. 
Therefore, local search is used to seek good feasible 
solution by changing widely the coefficient. 

I 1| 31 "I 5| 7| 8| 2| 8| 9| 12| 13| 10| 1l| 

i qi a ii ii ii jj t t 3 ' 3 ' 3 ' 4 ' 4 < 

Fig. 3 Representation of chromosomes (Upper: A series of 
work element numbers, Lower: A series of work center 
numbers). 



ol ol 'I D 



1 



H 2 I H 3| 3| 3| 4| 71 



I ol 11 ol ol 'I ol ol 'I ol ol 'I ol 'I 

Fig. 4 Representation of a series of work center numbers 
changed for crossover and mutation ( Upper: An original 
representation, Lower: A representation changed for 
crossover and mutation ) 
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3. Numerical Experiment 

3.1 Numerical conditions 

In order to evaluate the performance of the 
proposed algorithm, process design is performed on 
simple model by using the proposed algorithm. In this 
paper, process design for mix production is discussed. 
The model is constructed based on mix production 
including three types of products named as Product A, 
B, and C. Here, Product A is core product. Product B 
is a product based on core product and several parts 
are exchanged. Product C is a product constructed by 
including additional parts in the core product. Product 
A, B, and C are manufactured from 11, 11, and 14 work 
elements respectively. Each work element is determined 
as assembly of a part to semi-finished product. Table 
1 shows operation times of work elements. 'Av' and 
'#D' denote average of operation time and standard 
deviation of operation time respectively. Arrow to 
going to the left denotes that data in the column is the 
same as data in the left column. 

Fig. 5 shows procedure diagram of work elements 
of all products. The numbers in circles denote work 
element numbers in this figure. Character 'A', 'B', 
and 'C denotes Product A, B, and C. When work 



elements are similar to each other in mix production, 
the elements are combined to construct a single work 
element. Groups of work elements surrounded by 
polygons are regard as single work elements in the 
diagram. 

The amounts of products A, B, and C to produce in 
production plan are 100, 220, and 100 respectively. 
The period for production is 8400 min. (20 daysx7 
hours). Table 2 shows operation times in terms of 
groups of work elements regarded as single work 
elements, and these times are calculated from Eqs. 
(22) and (23). 

When multiple work elements are assigned to work 
center r, average operation time pi r and dispersion of 
operation time s r are calculated at the work center 
from Eqs. (22) and (23) considering the number of 
products. 



Mr 



Z^Z^J/Z^ 

^=zf^z%^]/z^ 



(22) 



(23) 



i V k Jl i 

Pi denotes the required amount of Product i. fi ik and 
s ik denote an average operation time and dispersion of 
operation time of work element k to produce Product i. 



Table 1 Operation time of work elements of products (No.: work element number, Av: average time, #D: standard 



deviation ) 


Product A 






Product B 






Product C 






No. 


Time 




No. 


Time 




No. 


Time 




Av 


#D 




Av 


#D 




Av 


#D 


1 


5 


0.3 


12 


8 


1.1 


1 






2 


5 


0.7 


13 


8 


1.6 


2 






3 


9 


2.1 


3 






3 






4 


9 


4.4 


4 






4 






5 


4 


0.4 


14 


7 


0.7 


5 






6 


7 


1.1 


6 






6 






7 


6 


2.4 


15 


9 


1.9 


7 






8 


4 


0.1 


8 






8 






9 


6 


2.2 


9 






9 






10 


5 


1.2 


10 






10 






11 


10 


0.2 


11 






11 
16 
17 
18 


3 
3 
4 


0.4 
1.3 
0.2 
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Product A, B, C 



(5')U 



Fig. 5 Precedence diagram of work elements of all 
products 

Table 2 Operation times of work elements combined by 
similarity of works ( No.: work element number, Av: 
average time, #D: standard deviation ) 



No. 


Av 


#D 


Work element number 


r 


6.57 


0.82 


1 (Product A,C), 12(Product B) 


2' 


6.57 


1.25 


2(Product A,C), 13(ProductP>) 


5' 


5.57 


0.58 


5(Product A,C), 14(Product B) 


7' 


7.57 


2.15 


7(Product A,C), 15(Product B) 



d ri k is 0 — 1 variable. When work element k to produce 
Product is assigned to work center r, <5„* takes 1. 
Otherwise, d„k takes 0. 

The minimum number of work centers is calculated 
from production plan by using Eqs. (1) and (2), and 
the number of work centers is determined to be five. 
Therefore, five work centers are assumed to be located 
in process line in this problem. As for genetic operator 
used in the proposed algorithm, both population and 
the maximum generations are determined to be 100. 

3.2 Evaluation of the Performance of the Proposed 
Algorithm 

Production lines designed by the proposed 
algorithm are compared with the lines designed by the 
different methods. Different methods and characteristics 
of designed lines are shown as follows: 

(Dl) A production line designed by using the 
proposed algorithm. Here, in terms of weighted 
coefficients to control values of criteria of infeasible 
solutions, the initial value of the weighted coefficient 
on left side, KL0, and the initial value of the 
coefficient on right side, KR0, are determined as 0.0 
and 20.0 respectively. The number of iteration to 



change weighted coefficient takes 15. 

(D2) A production line designed by the ranked 
positional weight technique for process design. 

(D3) A production line designed by resolving 
mathematical model for single objective functions 
presented in Section 2.1. The result is calculated by 
using free software "glpk.' 

The method D2 is used to design process for 
production line aimed at satisfying production plan. A 
purpose of this method is that work elements are 
assigned to work centers below target cycle time. 
Since this method tends to generate large difference of 
operation times among work centers, productivity of 
the line is ordinarily inferior to that of the line 
designed by optimal methods. In the process design 
obtained by Method D2, the maximum value of 
average operation time takes 19.42 and the maximum 
value of dispersion of operation time takes 23.77 
among work centers. 

On the other hand, the different lines designed by 
Method D3 give 15.99 as the maximum value of 
average operation time among work centers and 19.36 
as the maximum value of dispersion of operation time 
among work centers respectively. 

Fig. 6 shows solutions obtained by Method Dl. The 
distribution obtained by Method D 1 is one of popular 
and common results obtained by five trials. The figure 
shows distribution of all solutions in the final 
population and the solutions obtained by Method D2 
and Method D3. Filled circles and opened squares 
denote feasible solutions and infeasible solutions 
respectively. The solutions obtained by the proposed 
algorithm include the optimal solutions obtained for 
single objective functions by using Method D3. In the 
figure, several infeasible solutions are superior to 
feasible solutions and feasible solutions are generated 
at middle area between both axes. The proposed 
algorithm generates superior feasible solutions to the 
method of Method D2. 

In order to evaluate productivity of the lines, 
event-driven simulation is executed by using data of 
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production plannt methods 

work elements assigned at each work center by 
different methods, because different products flow in 
line for mix production and difference of operation 
times between successive work centers tends to be 
large. Fig. 7 shows comparison of productivity of the 
lines designed by different methods. Productivity of 
the line is evaluated by time to complete 420 products 
determined in production plan. It iterates the 
simulation 100 times. Solution 'A' and Solution 'B' 
generated by the proposed algorithm are evaluated 
using simulation in Fig.6. Fig.7 shows that the 
proposed algorithm is clearly superior to Method D2 
in the terms of productivity. These results show the 
proposed algorithm can generate Pareto and feasible 
solutions to provide high productivity for the problem. 
Since solution 'A' is smaller than solution 'B' in 
terms of maximum value of average operation times, 



average completion time of products obtained by 
solution 'A' is smaller than that obtained by the other 
solution. On the other hand, since solution 'B' is 
smaller than solution 'A' in terms of maximum value 
of dispersion of operation times, standard deviation of 
computational times obtained by solution 'B' is barely 
smaller than that obtained by the other solution. 

3.3 The Characteristics of the Proposed Algorithm 

In order to investigate the characteristics of the 
proposed algorithm in terms of the weighted 
coefficient, the proposed algorithm is compared with 
the single MOGA in the proposed algorithm under the 
conditions that the coefficient takes several static 
values between 1.0 and 50.0. Maximum generation 
takes 200 in order to take adequate time to seek Pareto 
solutions. Five trials are performed on each 
coefficient. 

When the weighted coefficient takes 1.0, the 
condition denotes that both chromosomes obtaining 
feasible solutions and infeasible solutions have the 
same probability to be selected in the process of 
crossover and selection for next generation in the 
single MOGA. When the weighted coefficient takes 
large number such as 20.0, 30.0, and 50.0, the 
condition denotes that chromosomes generating 
feasible solutions have higher probability to be 
selected for crossover and selection for the next 
generation than chromosomes generating infeasible 
solutions. 

Fig. 8 shows the ratio of feasible solutions to all 
solutions obtained by the single MOGA utilizing 
constant coefficients at the maximum generation. 
Each bar denotes average of the resultant ratios 
calculated from five trials and each error bar denotes 
standard deviation of the ratios. This figure shows that 
it is difficult to seek feasible solutions when the 
weighted coefficient takes less than 10.0. Therefore, it 
is difficult to seek high quality feasible and Pareto 
solutions under these conditions. When the weighted 
coefficient takes 1.0, no feasible solution is sought. 
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coefficient of 50.0 is larger than that obtained from the 
other coefficients. In this figure, some feasible 
solutions are located at areas far from Pareto front 
when the coefficient takes 20.0 or 50.0. We consider 
that the resultant distribution is obtained by easily 
selecting feasible solutions from bad feasible solutions 
when the coefficient takes large numbers. 

However, every coefficient generates similar 
solutions. This result denotes that the single MOGA is 
capable to seek high quality feasible and Pareto 
solutions when the coefficient takes large number. 

The weighted coefficient which generates adequate 
quality feasible solutions is not clear ordinarily. 
Therefore, control of the coefficient included in the 
proposed method is effective to seek feasible and 
Pareto solutions for problem generating many 
infeasible solutions such as process design. 
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Fig. 8 The ratio of feasible solutions to all solutions 
obtained by the single MOGA 
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Fig. 9 Distributions of all feasible solutions obtained by 
the single MOGA at the maximum generation when the 
weighted coefficient takes different numbers. 

On the other hand, when the weighted coefficient 
takes greater than 10.0, the ratio takes larger numbers. 
This result denotes that solutions are easily selected 
from bad feasible solutions when the coefficient takes 
large number. Therefore, we consider that local search 
to widely change the coefficient is useful to prevent 
from selecting solutions from bad feasible solutions. 

Fig. 9 shows distributions of all feasible solutions 
obtained by the single MOGA at the maximum 
generation when the coefficient takes 5.0, 20.0, and 
50.0. This figure shows the results in which the ratio 
of feasible solutions to all solutions takes 
approximately average number in five trials. The 
number of feasible solutions obtained from the 



4. Conclusions 

This study treats multi-objective problem for 
process design for mix production. Infeasible 
solutions are easily generated more than feasible 
solutions in the problem because precedence 
relationship of work elements does not easily satisfy 
by changing assignment of the elements slightly. 

New type of a multi-objective genetic algorithm is 
proposed and the proposed MOGA includes local 
search including control of the probability to select 
feasible solutions while both feasible solutions and 
infeasible solutions are stored. The proposed MOGA 
is performed on simple model of process design for 
mix production to evaluate the algorithm. Numerical 
experiments show that the proposed algorithm is 
effective to design production line which provides 
high productivity as well as high quality feasible and 
Pareto solutions. In addition, the results show the 
effectiveness of control to select infeasible solutions 
in the proposed algorithm. 

In future study, the proposed algorithm is applied to 
large scale and complex problem of process design. 
Furthermore, search process in the single MOGA is 
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developed as PSO and MOGA including specific local 
search to effectively generate high quality feasible. 
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1. Introduction 

Vehicle testing involves fabricating different vehicle 
parts that should be tested in different aspects. A 
vehicle has lots of constraints, for instance emissions, 
noise, space, handling and looks. This testing normally 
takes some months. 

A vehicle is defined by millions of parameters. 
Volvo can not do one test vehicle per defining 
parameter. Due to the high number of parameters and 
the long response time, the Lean Optimization method 
has been introduced. The Lean Optimization method 
combined with Linear Programming of the meta-model 
seems to be an efficient tool in optimizing vehicles. 

Lean Optimization method is a method introduced 
by Iana Siomina and Sven Ahlinder. It is based on the 
science of using experiments, Design of Experiments, 
DoE. In DoE, a planned and designed experimental 
series is performed and the steepest ascent is calculated. 
In DoE the number of experiments often is larger than 
the number of parameters. In Lean Optimization 
however, the number of parameters (162 in our test 
problem) is far more than the number of experiments 
for each step (14 in our case). This means that the 
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gradient is not calculated for each parameter in the 
problem one at a time, but is rather to be estimated in a 
sub space of the problem space. This method needs 
much less evaluations of the problem which is 
important if the evaluation time is days or years and the 
parameters are thousands or millions. 

Sven Ahlinder and Ivar Gustafsson concluded that 
the expected angle v between the estimated and true 
gradient depends on the number of rows, m, and 
columns, n. in the design matrix. E(cos(v))=square 
root(m/n) for large m, m<n. Assuming m/n=l/4 gives 
cos(v)=l/2 so the step taken in the approximate 
gradient direction gives half the gain compared to a full 
investigation, but allocates only 25% of the resources. 

For this problem, m=162 and n=14 which gives a 
ratio of gain=sqrt(14/162)=29%. This means that 
spending just 9% of the efforts we gain 29% of the 
"correct" steepest ascent. The normal, gradient, 
estimation of the ascent is calculated by 162 
experiments and will give 1 00% true ascent direction. 

In this article constraints on both in-data and 
out-data are considered trying to converge to an 
optimal solution of a constrained problem. The linear 
meta-model of the problem together with the linear 
constraints on in and out data is fed into a Linear 
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Programming solver and a good solution is gained 
using only a few evaluations of the probably expensive 
function. 

The problem is to optimize a physical object with 
several limitations. Since the method is not still 
established it could not be tested on a real physical 
object. The example studied here is a genuine Volvo 
calculation program for designing hardware for a 
heavy duty diesel engine. This calculation program 
behaves in many ways as a physical engine. That is 
why it is so valuable both for engine development and 
for this test. This engine calculation program is called 
"gas exchange calculation" and is very rapid (minutes) 
in response. Other calculation programs, 
Computational Fluid Dynamics for instance, have 
response times of weeks on a super computer. 

In this calculating program, emulating a Volvo D 1 3 
engine, there are some discrete variables, but the over 
100 variables considered are continuous and limited 
with a lower and upper value. Those are for instance 
bore and compression ratio of the cylinder, diameters 
and lengths of all tubing, timing and other injection 
parameters and parameters for the turbo and EGR 
systems. The out-data considered is brake mean 
efficient pressure, BMEP, the specific power of the 
engine, which is maximized with constraints on 
emissions and forces of the engine. All evaluations are 
made at maximum injected fuel amount for the 
considered engine. 

The optimization of the BMEP could have been done 
using a normal gradient method since the response time 
for the program is only a few minutes. The Lean 
Optimization method, however, is aimed for 
computational fluid dynamics programs, CFD, with 
response times of days or weeks. There are even 
hand-made physical engines with response times of 
months or years where this can be applied. The task is 
to do this optimization using as few evaluations as 
possible, which becomes important when performing 
experiments where the response times are weeks or 
months. 



2. Method 

The Lean Optimization method has showed a good 
potential on synthetically problems [1]. It also has a 
good theoretical potential [2]. It has though neither 
been applied on a real construction problem, nor been 
applied to problems with constraints in the responses. 

The Lean Optimization method is related to Design 
of Experiments, DoE, [3] since in both cases an 
experimental series is run on the problem. The idea of 
DoE is convenient since it is defining a batch of 
experiments that could be prepared simultaneously. 
This is convenient both in calculations and using 
physical items when all parts needed could be produced 
before the series of experiments begins. 

In gradient method all parameters are changed one at 
a time in a small interval delta. Suppose the response of 
the changes could be written f(X), then X could be the 
matrix, table, where each row is one experiment and 
each column is one variable. X will consist of a row of 
zeros for the original point and under this, a diagonal of 
deltas. Normally each row will be evaluated a separate 
experiment but it is equally true to consider X as an 
experimental series with f(X)=Y, where each column 
of Y is the vector of all results, variable by variable of 
the experiments in X. Regarding the gradient as a 
system of equations XB=Y, where each row of B is the 
coefficients of the gradient, the ascent direction. 

In Design of Experiments, DoE, X is not a diagonal 
matrix, but a balanced orthogonal matrix consisting of 
"-1" and "+1" [3]. The cause of this is to reduce noise 
influence in the experimental series. The DoE method 
gives the system "a solid kick" [3] which means that a 
reasonably large interval in the in-parameters is 
spanned by the experimental series. 

DoE has become popular in numerical calculating 
since the system could not only be approximated by 
linear models, but also with interactions, xl*x2, and 
quadratic terms. Having a system with interaction and 
quadratic terms, this may lead to an optimum in fewer 
experimental series but the total numbers of 
evaluations, experiments, is using just linear terms 
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assumed fewer. 

A Lean Optimization design could be made in the 
following way: Take a full factorial two level design 
[3] , all against all in two levels. This could be a design 
matrix consisting of 10 columns and 1024 rows. 
Transpose this table, matrix, so the observations are 
used as in-parameters and the parameters are used as 
observations. This gives a short and fat experimental 
design X. The first half of the columns will be identical 
to the last part of the columns except for the sign. 
Therefore [2] recommends to cut off the last half of the 
transposed full factorial design and then cut off a 
suitable number of columns more to get the right 
number of columns of variables. The matrix X 
eventually consists of 10 rows (experiments) and 
maximum 512 columns (variables). This design is 
called Ahlinder-Gustafsson-design AG-design [4] and 
has more or less orthogonal rows (experiments) which 
will be shown useful. 

In this work, a matrix of 14 rows and 162 columns 
were used. The reason of the large number of 
experiments were due to a wish that every columns 
should not only be independent of the others but also be 
balanced [3], which means there are equally many "-1" 
as "+1" in each column. This is used for minimizing the 
influence of noise and may not be applicable in a 
calculation program. 

In the Lean Optimization method, the experimental 
design X is a super saturated design [5] which from an 
evaluation point of view will give an under-determined 
system of equations. This system of equations is solved 
using minimum norm solution of Moore -Penrose 
inverse [2] invX=X'(XX')-l. In Matlab, B is easily 
calculated as B=pinv(X)*Y.B consists of the columns 
of ascent directions that are not necessarily the steepest 
ascent but has an expected angle v to it so 
E(cos(v))=square root (m/n) for large m, m<n, where m 
is the number of rows or experiments and n is the 
number of columns or parameters studied [2]. This 
implies that if m/n=l/4, cos(v)=l/2 and the step in 
ascent direction taken gives half the gain compared to a 



full investigation, but allocates only 25% of the 
resources. 

There are several articles written about "super 
saturated design" [5]. In these articles, generally the 
aim is to make the columns, not the rows, of the design 
matrix as independent as possible. This is due to that 
the system of equations is solved using variable 
selection [3] not minimum norm, as is done here. 

To handle constraints of emissions and forces, a new 
method was introduced. It uses linear programming, 
LP [6]. Linear programming solves linear problems 
with linear constraints. If the experimental design is 
solved for all out-data, Brake Mean Effective Pressure, 
BMEP, emissions and forces, then we have a linear 
(not perfect) meta-model of our problem. All responses 
in an evaluation point, the vector e, consisting all 
in-parameters will form a linear function g(e)=e*B, 
where e*B is a set of linear functions and each element 
in g is a response. The first element will be the BMEP 
that should be minimized fulfilling constraints on the 
other elements in g, the emissions and forces. 

BMEP is the linear target function, emissions and 
forces are linear constraints on out-data and the 
investigated area is a set of linear constraints of the 
in-data since the in-data has been given a solid kick. 
Somewhere in this linear constrained space, there is an 
optimum of the linear function of BMEP fulfilling the 
constraints and the LP solver will find it. 

This reached optimum for one experimental design 
is followed by another step of experimental design 
centered around the first optimum. The optimization 
will continue in several steps each reaching an 
optimum for BMEP in the same way as the gradient 
search does. 

3. Results and Discussion 

The Lean Optimization method combined with an 
LP solver works well on the engine design program. 
This means that it will work well on equally complex 
problem which could be CFD calculation or even 
physical engines. The LP solver finds the optimal 



Constraints in Lean Optimization Method Tested on a Real Engineering Problem 



223 



solution of the linear problem which is the result of the 
Lean Optimization approach. 

The result of the optimization series can be seen in 
Fig. 1. Each dot is a step in the Lean Optimization 
method. The method has converged somewhere 
between 100 and 150 evaluations. This means that the 
Lean Optimization method in this (real) case has 
converged before the gradient method would be able to 
take its first step, since it has to do 162+1 evaluations 
for each step for the 162 variables. 

According to Fig. 1, the engine could be downsized 
8% to increase BMEP 8% with constant emissions and 
other constraints. The piping, injection and turbo 
system was well optimized already thanks to hard 
working engineers. The industrial solution was to 
increase the power of the engine 8% instead of 
remaking the whole engine design. 

The Lean Optimization method projects the problem 
on a linear subspace, which rank is dependant of the 
number of rows or experiments in the design. Both 
target function and constrain functions are projected on 
the investigated subspace. This gives linear models of 
both target function and constrain functions. Probably, 
LP solver is useful even for full rank gradient or DoE 
methods. In DoE, there are also non-linear models 
which do not give linear problems 

In the relative low dimensional model with 162 
in-parameters, a standard LP solver (Matlab) will work 
well. The LP solver will find an optimum in a corner of 
the convex envelop consisting of constraints of in-data 
and responses where BMEP has its maximum. This 
point is chosen as starting point of the next 
investigation. 

It seems by the testing that using a new subspace for 
the new step in the ascent direction is more efficient 
than using the same subspace that was used in the 
previous step. This is made by recombine the columns 
in the lean design matrix. The new subspace could be 
randomly or analytically recombined. The difference is 
not measurable, possibly due to the low rank of the 
subspace. For the 162 variable problem only 14 
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Fig. 1 Optimization result showing that Lean Optimization 
converges in less number of evaluations than the number of 
variables (162). 

experiments were conducted each time. This makes the 
random subspace quite independent of the previous 
subspace. However, using the same subspace twice 
seems as a poor approach. This is possibly due to that 
the potential of this subspace is exhausted by the step 
taken. 

When the optimum points backwards, two 
consecutive directions scalar product is negative, the 
search area was shrunken. This is generally done in 
DoE optimization and works well here also. 

The result was that BMEP could be increased 8% 
just shrinking the displacement of the engine. The other 
over 1 00 variables was more or less unchanged. This is 
natural since those variables have been optimized by 
engineers in generations. The remarkable thing is that 
this conclusion was drawn without applying any 
engineering judgment, just mathematics. The result is 
used to increase the power 8% of the existing engine. 

4. Conclusions 

Lean Optimization combined with linear 
programming seems to be a powerful tool in 
constrained optimization of complex functions such as 
engine hardware. It is demanding to calculate each step, 
but when the evaluation time of the problem is days or 
years, the goal has to minimize the number of 
evaluations, not the computing time. 
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Physical models and physical objects could be 
non-linear and not well approximated with a linear 
model. In the engine design case, the physical based 
model is reasonably linear so it is useful with a linear 
approximation. 

Shrinking the search area when having a negative 
consecutive directions scalar product, is a useful 
method also in Lean Optimization method. 

It seems vital to permutate the columns of the lean 
design matrix to search new dimensions of the 
parameter space for each step. 

In this example, the Lean Optimization method 
converges before the gradient method would have 
taken its first step.) 
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(b 0 + 6 1 t)x(t — 1) + (ci 0 + d^xit + 1) + /(t), t £ R, (*) where f(t) — En=o fnt n - This equation is investigated with the use of 
the method of polynomial quasisolutions based on the representation of an unknown function in the form of polynomial x(t) — 
2™ =0 x n t n . As a result of substitution of this function into equation (*), there appears a residual A(t) = 0(t w ), for which an exact 
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1. Introduction 

Investigation of any process that occurs in the 
modern world starts as a rule with creation of its 
physical model. Further this model is assigned a 
mathematical model which is studied by various 
mathematical methods. In many cases the 
mathematical models represent the systems of 
differential equations. However, the structure of these 
equations can be so complex that it is not possible to 
apply directly the methods available for studying the 
obtained equations. Then, we transfer from the initial 
system to the modified system which allows the 
application of available methods to study it. For 
example, the study of nonlinear equations employs 
linear equations, variable parameters of equations are 
replaced by constant parameters, etc. The results 
obtained on the basis of modified mathematical model 
require additional studies from the view point of their 
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correspondence to solutions of physical model. In 
applications a great role is played by ordinary linear 
differential equations for which a great number of 
research methods have been developed. But, unlike 
rather satisfactory state of the theory of ordinary 
differential equations detailed research into the 
processes that occur in the world around makes us 
investigate more complex equations and take into 
consideration the fact that variables depend not only 
on the states at a current time instant but on their 
values at the previous and subsequent time instants as 
well. Besides, many problems lose sense if the 
relationship between variables and their values at 
different time instants are not considered. In these 
cases we use functional differential equations as 
mathematical models. The equations that take into 
account prehistory of the processes first appeared in 
the studies conducted by I. Bernoulli and L. Euler. At 
the beginning of the 1940s- 1950s the interest of 
mathematicians in the research on such differential 
equations increased greatly owing to the problems of 
the control system theory and also because these 
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equations have very interesting properties. 
Nevertheless, this area intensively developed only in a 
very limited number of aspects and by a rather limited 
number of mathematicians. An impetus to the study of 
delay differential equations by the mathematical 
methods was given in the monograph by A.D.Myshkis 
[1], that was published in 1951. Here, it is necessary 
to note successive fundamental studies performed by 
E.Pinni [2] R.Bellman and K.L.Kuk [3], N.V. Azbelev, 
V.P. Maximova and L.F.Rakhmatulina [4]. Currently 
owing to the efforts of mathematicians from different 
countries the theory of differential difference 
equations is thoroughly developed in different areas, 
namely in the theory of automatic control, the theory 
of self-oscillating systems, in the studies related to 
fuel combustion in the rocket engine, in long-term 
forecasting in the economics, in some biophysical 
problems and in many other fields of science and 
technology, whose number is persistently growing. 
The abundant applications encourage rapid 
development of the theory of functional differential 
equations (FDE) and currently this theory belongs to 
the most developing areas in mathematical analysis. 
Nevertheless, even in this seemingly well studied area 
still there are many problems solutions to which 
present a great interest in terms of both science and 
application. 

This paper focuses on the linear differential 
difference equations (LDDEs) of mixed type. 

x(t) = a(t)x(0 + b(t)x(t - 1) + d{t)x{t + 1) + 
f(t),tER. (1) 

Interest in the equations that contain simultaneously 
delay and advanced arguments is caused by a number 
of problems in which these equations appear as a 
necessary object for study [3], [5], [6]. When 
coefficients of the equation are constant or have some 
special representations, application of Euler's method 
allows one to obtain analytical solutions, which are 
generated by roots of the characteristic 
quasi-polynomial. The authors have no information on 
the solvability conditions for LDDEs in the class of 
analytical functions in the case when parameters of the 



equation are variable. When coefficients of the 
equation (1) are represented by polynomials, a method 
of polynomial quasi-solutions (PQ-solutions) has been 
proposed in [7]-[10]. It implies that a formal solution 

n _ Q X„t" is 

introduced into the investigation. Hence the term 
"polynomial quasi-solution" (PQ-solution) is 
understood in the sense that after its substitution into 
the initial-value problem there appears a residual 

A(t) = 0(t") , for which an exact analytical formula 

has been obtained. On the other hand, this allows one 

to find unknown coefficients x n and, consequently, 

the PQ-solution itself. The paper is devoted to finding 
out the conditions for the existence of PQ-solutions to 
the initial-value problem with the initial point for 
some LDDEs of mixed type, as well as to the 
techniques of their obtaining. 

2. Problem Statement 

Consider the following initial problem for scalar 
linear differential difference equation of mixed type. 
x(t) = a(t)x(t) + b(t)x(t - 1) + d(t)x(t + 1) + 
f(t),teR,x(0) = x o , (2) 

where 

a(t) = a 0 + a-tt, b(t) = b 0 + b t t, d(t) = d 0 + 

d l t,f{t)=YJn=ofnt n . (3) 
Introduce the polynomial 

x N (t) = ^ =0 x n t n ,teR. (4) 

In this case 

x N(t) = £n=o nx nt n 1 , (5) 

X w (t — 1) = £n=o x n(.t — 1)" = £n=o ' < (6) 

x N (t + 1) = E^o x n {t + l) n = Z5U x nt n . (7) 
Here 

j _yN-n r i 

t _ (n + Q! _ _ ^ 

C-n+i ~ ^T^j > n — 1,N — 1,X N — X N — X N . 

We analyze the dimensions of polynomials 
obtained by substituting (3)-(7) into (2). The 
derivative has the dimension N — 1 , the terms 
a(t)x(t) , b(t)x(t - 1) , d(t)x(t + 1) have the 
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dimension N + 1 , and /(t) - the dimension F. On 
the other hand for the last coefficient x N in (4) to be 
determined by the last coefficient f F in (3), while 
applying the method of undetermined coefficients, it 
is necessary that N = F + 1 in (4). Define the 
function f(t) as 

/(0= In=o/nt n + A(t), (9) 
where fi=fi,i = 0,F = N — l , are known 
coefficients ( at the same time the cases where some 
or all fi = 0 are not excluded), and residual 



A(t) = f N t N + f N+1 t N+1 , f N and f N+1 are 
unknown coefficients. Based on the introduced 
notations we consider the initial problem 

x(t) = a(t)x(0 + b(t)x(t - 1) + d(f)x(t + 1) + 
f(t),teR,x(0) = x o . (10) 

Definition 1. Call problem (10) coordinated in 
dimension of polynomials with respect to problem (2). 

By substituting (2)-(7) and (9) into (2), by the 
method of undetermined coefficients we obtain 



nx n = \ ll =0 O^n-l-i + Mn- 



fo.n = 1; 

+ df*n-i-f) + fn-i, 2<n<N; 



0 = I;=o (.a-i^N-i + biX N -i + di^N-i) +f N ,n = N + l; 



(11) 



(12) 



0 = (a 1 + b ± 

Remark 1. Since the degree of polynomial x(t) is 
equal to F + 1, this makes it possible to select the 
degree of polynomial f(t) in (3) depending on the 
required degree of polynomial x(t), by adding a 
corresponding number of zero members to f(t). 

Definition 2. If there exists a polynomial of degree 
N = F + 1 

x N (t) = Y,n= 0 x n t n ,t£R, (14) 
that identically satisfies problem (10), then this 
polynomial will be called a polynomial quasi-solution 
(PQ-solution) to problem (2). 

Thus, the objective is to identify the conditions for 
the existence and methods for finding unknown 
coefficients f N and f N+1 , that make it possible to 
represent a solution to the initial problem (10) in the 
form of polynomial (14). 

3. Theorem of Existence of Polynomial 
Quasi-Solutions 

According to definition 2, the issues about the 
existence of polynomial quasi-solutions consist in 
establishing the conditions, under which we can 
calculate the coefficients x n , n = 0, N of polynomial 
quasi-solution (14) and coefficients f N and f N+1 of 
residual A w (t). To this end, let us consider relations 



+ dJx N + f N+1 ,n = N + 2. (13) 

(1 1) (13). Assuming that n = 1,2, (F + 2) and 
taking into account (8), we rewrite these formulas in 
the following form 
(-b 0 + d 0 - l)x 1 + (b 0 + d 0 )x 2 + (-b 0 + d 0 )x 3 
+ - + «-l) N b 0 + d 0 )x N = 
-(a 0 + b 0 + d 0 )x 0 - f 0 ; 
(a 0 + b 0 -b 1 + d 0 + di)*! + (C%(-b 0 + d 0 ) + b 1 

+ d 1 - 2)x 2 + 
(Ci(b 0 + d 0 )-b 1 + d 1 )x 3 + - 

+ «-l) N - 1 (Ct 1 b 0 -b 1 ) + d 0 
+ d ± )x n = 
+ b ± + djxg - /i; 
(% + b 1 + d 1 )x 1 + (a 0 + b 0 + d 0 + C^—b-L 

+ d 1 ))x 2 + (Cl{-b 0 + d 0 ) + 
CiQh + d ± ) - 3)x 3 + - + ((-l)"" 2 ^" 2 ^ 

- ct'h) + 

C$ 2 d 0 + C$ 1 d 1 )x N = —fzl 

(% + b 1 + d^Xiv-2 + (a 0 + b 0 + d 0 + C^.^-^ 
+ d ± ))x N _ ± + 
(Ch(-b 0 + d 0 ) + C 2 ^ + d ± ) - N)x N = -f F ; 
(a ± + b ± + d ± )x N _ ± + (a 0 + b 0 + d 0 + C^-b-^ 
+ d{))x N = —f N ; 
(a ± + b ± + d ± )x N = -f N+1 . (15) 
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These relations represent a system of linear following matrix form: 



algebraic equations of relatively unknown x n ,n = 

0, N and coefficients f N and f N+1 which has the 

/-b 0 + d 0 -l 
a 0 + b 0 - b 1 + d 0 + d ± 
a 1 + b 1 + d x 
0 
0 



M„ = 



0 
\° 



(-l) w 6 0 + d 0 

{ _ 1)N -2 {C N-2 bQ 



M N x^ — g* N , 

where 
b 0 + d 0 

C 2 i (-b 0 + d 0 ) + b 1 + d 1 -2 
a 0 + b 0 + d 0 + C\{-b x + d ± ) 
a 1 + b 1 + d x 
0 



0 
0 
0 

■b ± + C 
1 



(16) 



r N- 
L N 
r N- 

■N- 
N 



1 d 0 + d ± ) 

2 dn + 



N 
N 

i) + Cff 

N-3 



r N 
^1 



%) + cfl- 3 d 0 + C^~ 2 d 1 
%) + cfl- 4 d 0 + C^~ 3 d 1 



-Ch(.b 0 -d 0 ) + Cfi(b 1 + 2d 1 ) 
a 0 + b 0 + d 0 - C^(&! - dj) 



0 
0 
0 
0 
0 

0 
1 
0 



0 
0 
0 



(17) 



Z N — ( x l> x 2> — > X N> fN> fN + lV 



(18) 



5w = (-(«o + b 0 + d 0 )x 0 - f 0 , -(% + b ± + 

d 1 )x 0 -f ll -f 2 -f F ,0,0) T . (19) 

The following theorem is true. 
Theorem 1. Let matrix M N of linear system (16) 
which is determined by virtue of equation (2), be 
nonsingular. Then for any x 0 6 R the initial problem 
(2)-(3) will have a polynomial quasi-solution in the 
form of polynomial (14) of degree N = F + 1 with 
residual A(t) = f N t N + f N+1 t N+1 . 

Indeed, in this case there exists the inverse matrix 
N and from (16) we obtain 



teR,x(Q) = l. (20) 
Taking into account (14), we determine the 
polynomial quasi-solution as follows 



x n(0 — ^ x n t n . 



of matrix M 



**n = Mn^n- 



Hence, coefficients x 1 ,x 2 ,—,x N of polynomial 
quasi-solution (14) and coefficients f N and f N+1 of 
residual A(t) are calculated uniquely, which proves 
the theorem. 

4. Numerical Experiment 

Consider the initial problem with an initial point for 
the following LDDE: 

x(t) = (-2 + O.St)x(t) + (1 + O.St)x(t - 1) - (2 
+ 0.5t)x(t + 1), 



n=0 

Then, according to definition 1 , we write the initial 
problem coordinated in polynomial dimension 

x(t) = (-2 + 0.5t)*(0 + (1 + O.St)x(t - 1) - (2 
+ O.St)x(t + 1) +A N (t), 
t G R,x(0~) = 1, (21) 

where A w (t) = f N t N + f N+1 t N+1 . 

Consider an algorithm for finding polynomial 
quasi-solutions with N = 4. In this case linear system 
(16) will be written as follows: 
M 4 x* 4 = g\. 

Here in accordance with (17) - (19) 

/-4 -1 -3 -1 
/_4 _8 -4 -12 
0.5 -0.5 -12 -10 



(22) 



Ma = 



0.5 

0 
0 



-6 -16 
0.5 -7 
0 0.5 



0 
0 
0 
0 
1 
0 



\ 



V 



Polynomial Quasisolutions Method for Some Linear Differential Difference Equations of Mixed Type 229 



X4 — (Xi, x 2 , x 3 , x 4 , f 4 , fs) T , 
g* 4 = (3,-0.5,0,0,0,0) r . 
Since detM 4 = 5120 ^ 0, we find the inverse 
matrix 



M4 1 



n 



-0.2241 

1023 
-0.0795 
0.0330 
y 0.2708 
\-0.0165 
Then 



x\ = M 4 1 g\ 



-0.0170 

-0.1072 

0.0680 

-0.0289 

-0.2360 

0.0144 



0.0703 

-0.0391 

-0.0918 

0.0332 

0.2783 

-0.0166 



-0.0172 

0.0984 

0.0113 

-0.0637 

-0.4514 

0.0318 



0 \ 



V 



= (-0.6642,0.3606, -0.2725,0.1135,0.9304, -0.0567) T . 
Hence, we obtain 
x 4 (t) = 1 - 0.6642t + 0.3606t 2 - 0.2725t 3 
+ 0.1135t 4 , 
A 4 (t) = 0.9304t 4 - 0.0567t 5 . 
Let us present the results of calculations of 
polynomial quasi-solutions for N = 5 — 9. 

x 5 (t) = 1 - 0.7134£ + 0.3378t 2 - 0.1560t 3 
+ 0.1309t 4 - 0.0497t 5 , 



A 5 (t) = -0.4631t 5 + 0.0249t 6 , 
x 6 (t) = 1 - 0.7264t + 0.4027t 2 - 0.1310t 3 
+ 0.0572t 4 - 0.0609t 5 
+ 0.0215t 6 , 
A 6 (t) = 0.2241t 6 -0.0108t 7 , 
x 7 (t) = 1 - 0.7077t + 0.4206t 2 - 0.1911t 3 
+ 0.0336t 4 -0.0167t 5 
+ 0.02279t 6 - 0.0093t 7 , 
A 7 (t) = -0.1065t 7 - 0.0046t 8 , 
x B (t) = 1 - 0.7000t + 0.389t 2 - 0.2175t 3 

+ 0.0836t 4 + 0.0030t 5 + 0.0023t 6 
- 0.0128t 7 + 0.0040t 8 , 
A 8 (t) = 0.0509t 8 - 0.0020t 9 , 
x 9 (t) = 1 - 0.7100t + 0.3650t 2 - 0.1777t 3 

+ 0.1153t 4 - 0.0353t 5 - 0.0130t 6 
+0.0018t 7 + 0.0061t 8 - 0.0018t 9 , 
A 9 (t) = -0.0246t 9 + 0.0009t 10 . 
The plots of polynomial quasi-solutions are given in 
Fig. 1. 
40 T 



x 4 (0 
x 5 (t) 
x 6 (t) 
x 7 (0 
x s (t) 
x 9 (t) 




Fig. 1 Polynomial quasi-solutions 
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Let us introduce the following definition. 
Definition 3. s -attractability of polynomial 
quasisolutions at some interval [t 0 , t ± ] is 
understood as the property of reciprocal attraction of 
a sequence of polynomial quasisolutions generated 
by an increase in degree N of the polynomial of 
polynomial quasisolutions, i.e. there exists N* such 
that for N > N t and for a given s 
\x N+i (t) -x N+ ;-i(t)| <£,i = 1,2 ...k,Vt 6 [t 0 , til- 
According to this definition, for e = 0.01 
polynomial quasi-solutions x n (f) of problem (20) at 
interval (—1,1) possess the property of 
£-attractability. 
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interviewed all over the country with the objective of studying relationship between different variables such as disability, level of 
education, gender, income in minimum wages among others. 
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1. Introduction 

It is currently regarded as a fact that has always 
existed throughout history people with disabilities [1, 
2]. Gradually, societies perceive that beyond the 
charity and assistance, such persons should be 
included in programs and policies that could enhance 
their productive potential [3, 4]. In fact, the very 
people with disabilities were showing signs that they 
could and wanted to study work and be fully included 
in society [5]. 

Thanks to the mobilize these people, it was possible 
to secure nowadays, a set of laws that supports 
populations, not only with regard to the labor market, 
but also as a human rights such as education, health, 
leisure, and finally, the right to exercise full 
citizenship [6]. 

According to the World Health Organization (WHO) 
estimates that more than a billion people live with 
some form of disability, something close to 15% of 
the world population (based on 2010 estimates). This 
is lower than the estimates coming from the IBGE, 
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which date back to 2010 and suggest that about 23.9% 
of Brazil's population, are people with at least one 
disability. 

Even with these advances, whether in historical 
terms, as the paradigm shift that allow people with 
disabilities to be citizens, as in the legal field with the 
existence of laws and decrees that established these 
rights, however, the participation of these people also 
is very restricted. According to the IBGE Census 2010, 
it is estimated that a little less than 5% had completed 
higher education, which indicates that a high number 
of people with disabilities continue exerting informal 
activities, poor and discontinuous or simply have no 
occupation, living based on retirement, pension and or 
family support. 

It is believed that low working conditions of people 
with disabilities are due to situations such as difficult 
access to education, inadequate infrastructure, lack of 
information and bias on the part of schools and 
companies that make these people with a lower level 
of schooling impeding the entry of them in the formal 
labor market. 

According to WHO [7] in regards to disabled 
persons notices as: 

(1) Population of people with disabilities in greater 
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growth driven by aging; 

(2) Diverse experiences resulting from the 
interaction between health conditions, personal and 
environmental factors; 

(3) Most vulnerable populations indicating higher 
prevalence in lower income countries; 

(4) Presentation levels of poorer health; 

(5) Lower educational achievement with lower rates 
of stay and lower school performance; 

(6) Lower economic participation for people with 
disabilities more likely to be unemployed and earn 
less when employed; 

(7) Higher rates of poverty, suffering high rates of 
deprivation with food insecurity, poor housing, lack of 
access to clean water and sanitation and inadequate 
access to health services, and finally; 

(8) Greater reliance and participation restrictions 
due to reliance on institutional solutions, lack of 
community life and inadequate services. 

The reverting this picture is no simple task. 
According to experts [5], this requires actions by at 
least eleven key aspects: 

(1) Expanding public knowledge about people with 
disabilities and their integration in education and 
work; 

(2) Issues of law (not only regarding the so-called 
"Quotas Law", but also the labor laws / social 
security); 

(3) The strengthening of educational inclusion and 
possibilities of professional qualification and decrease 
the deficit amount of people with disabilities to be 
hired to meet the " Quotas Law " [6,8]; 

(4) More incentives to entrepreneurship triad in 
encouraging more free enterprise, innovation in order 
to encourage more creative solutions to a potential 
consumer market and operation of small business 
ideas to business, opening and closing company, niche 
markets in order to offer a product with a high 
probability of being well received and as forms 
capture the potential customers and have good output 
from potential consumers. For this, we believe that 



these notions points should be part of school curricula 
since the fundamental level so that these ideas are 
more mature in the minds of students in subsequent 
levels and complemented by the government and 
private enterprise with facilities for opening small 
business more cost subsidized financing more 
facilitated with lower interest rates offered by the 
government and financial system and longer terms of 
payment and existence more easy media entities 
linked to small and medium enterprises as SEBRAE 
and innovation as the USP, UNICAMP and IPEN; 

(5) Accessibility as a concept summary of an 
inclusive society that much more than adequate 
infrastructure, creating conditions so that people with 
disabilities can use products, services and information 
as any other citizen; 

(6) The consolidation of new paradigms and ways 
of thinking about disability issues in society in general, 
but especially between schools and employers and 
people with disabilities themselves; 

(7) It should be considered mandatory for health 
professionals and care social notify the Ministry of 
Health and met people who have been diagnosed with 
problems related to their disability with CIF's; 

(8) establishment of a quota of 5% of seats for 
admission of persons with disabilities in schools in 
different educational levels in order to stimulate a 
greater number of people with disabilities to acquire 
or improve their level of education; 

(9) Implement methodologies for collecting data on 
people with disabilities, standardizes them, let them 
internationally comparable so that they can establish a 
benchmark and monitor progress and monitor policies 
related; 

(10) strengthening the research on disability so 
there is increasing public understanding of the supply 
of information for policymaking dedicated to 
disability, and for the efficient allocation of resources, 
and ultimately, and, by end; 

(11) More than Quotas Law, to ensure that the 
disabled person is hired for functions that correspond 
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to your true potential or ability [5,9]. 

Besides all this, it is necessary for economic and 
social conditions of the country evolve positively. The 
accelerated economic growth, better distribution of 
income, with quality public services and effective 
social programs, among other things, are beneficial to 
everyone, including, of course, for those with a 
disability. For more incredibly that it seems there isn't 
a world "specific" people with disabilities. 

However, to better assess the need and the 
fulfillment of this juridical structure, it is necessary to 
better describe this group of people to find answers to 
questions like how many? Where do they live? How 
do they live? What are the implications that disability 
involves access to all of these people human services 
autonomously and full? In short, as the deficiency 
may influence quality of life. 

In terms of statistics, shows that there are few 
studies in formal terms, among which stands out the 
data obtained from the census, which allows questions 
as: How people with disabilities are distributed 
throughout the country? How to assess the access of 
people with disabilities in terms of the various 
services mentioned above? How is the evolution of 
disabled people by comparing them with those 
without disabilities? Different deficiencies are 
homogeneous. It is possible to form homogeneous 
groups. What are the variables that contribute most to 
the problems of deficiencies? There is a national 
registry of people with disabilities? Responding to 
these and other questions in statistical terms may 
possibly contribute to better support these people in 
order to be better assisted and resources are better 
managed and optimized by the actions of public 
policies in this area. 

The Census is the most complex statistical 
operation undertaken by a country, when they 
investigated the characteristics of the population and 
households in the country and is the most important 
reference source for the knowledge of the living 
conditions of the population in all locations. 



In Census 2010, we used two types of questionnaire: 
Basic Questionnaire that was administered to all 
households and Complete Questionnaire or Full 
Sample that was applied to all households selected for 
the sample, which besides research contained in the 
Basic Questionnaire, includes other characteristics of 
the home and search important information social, 
economic and demographic characteristics of its 
residents. 

Questions relating for people with disabilities were 
applied in Sample Questionnaire, which sought to 
identify the deficiencies visual, listen and motor 
activity, according their degree of severity, through 
the perception of the people about these difficulties 
and for those who have declared intellectual or 
intellectual disability. 

In terms of statistics been applied, firstly, the 
homogeneity test to check if the degree of severity of 
the types of visual impairment, of hearing and 
locomotion behaves similarly to the homogeneous or 
different levels of variables such as age, educational 
level and income. 

Then using the technique exhaustive AID groups 
were formed for the study population according to 
each incidence of disability or not considering their 
degree of severity, there is at least one disability and 
the amount of defects considered (0-4). This method 
been applied to the responses of those interviewed. In 
each case, for each disability groups have formed for 
the sample population selected for Brazil at large, by 
region and state. 

The exhaustive AID algorithm investigates all 
variables, all groups were scrutinized and it selects the 
one that presents the greatest association with the 
dependent variable. This is done for all independent 
variables and can be used in situations where the aim 
is to divide the population into different segments for 
a given criteria. In this case, the criterion chosen was 
the effect of various deficiencies studied. 

The main objective of this work is to study the 
relationship between deficiencies of the population 
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and their ways of working and studying. 

The variables used in this study were obtained 
directly from the Database of Sample of respondents 
who answered the complete questionnaire and were 
divided into full blocks considering the following 
topics: identity, family, work, instruction and 
disability. 

2. Motivation 

To be able to include people with disabilities, it is 
necessary, first of all, more accuracy estimate what 
would be the amount of people who find themselves in 
these conditions for each of the different disabilities, 
how they live and where they live, and an alternative in 
this case, was to consider the database obtained from 
census 2010 for the Sample Questionaire, and, 
according to this same census estimate, it is believed 
that there are 45,606,048 persons in Brazil with at least 
one permanent disability, representing approximately 
23,9% of all the population represented, in the same 
census, by 190,755,799 habitants. 

2. 1 Homogeneity Test 

This test consists in verifying whether a random 
variable behaves homogeneous in several 
subpopulations and fixing the sample size in each of 
these subsets and then selects one sample each. 

To calculate the expected values, assuming that 
there is homogeneity between subpopulations, is used 
for each cell (i, j), according to expression (1) below: 

total in the column /' 

e i j = n i x (1) 

overall total 

The total line «, indicates the sample size of the 
subpopulation i, whereas the quotient, total in the 
column j divided by the overall total, representing the 
proportion of occurrences of the variable value 
corresponding to the column j. If there is homogeneity 
of the behavior of the corresponding variable, it is 
expected that this ratio is the same in all 
subpopulations. 

The next step is to calculate Q 2 that is the difference 



between the observed and expected values using the 
expression (2) below: 

. On - e, 7 



(2) 

i=\j=\ "y 

where r and s represent the number of row and column 
respectively. 

To interpret the expression of Q 2 , we note that the 
term oy - ey indicates the difference between the 
observed value and expected value in row i and 
column j, if there was homogeneity [10]. 

For a large number of observations, the expression 
of g 2 is chi-square with (r - 1) (s - 1) degrees of 
freedom. The critical region contains large value of Q 2 , 
i.e., BC =Q 2 :Q 2 >q c , with q c , being determined by 
the significance level of the test, in other 



words, a = p[q 2 > q c /h 0 is true) 



In the case of this study, homogeneity tests were 
made among variables that represent the different 
disabilities and variables such as educational level, 
gender and income level. 

2.2. Technical exhaustive AID 

It is a technique used in situations in which the goal 
is to divide the population into differentiated segments 
in relation to a given criterion. The exhaustive AID is 
based on chi-square test on a contingency table 
between the categories of the dependent variable and 
the categories of the independent variables [11,12]. It 
is a statistical method for efficient segmentation. In 
this particular case, the criterion chosen was the 
incidence of defects in each of the study [13,14]. 

The final objective is to build a classification tree 
containing only the most important variables for the 
classification, along with their most significant 
categories for the response variable. The steps of the 
algorithm exhaustive AID are following: i) definition 
of the research problem; ii) characterize the sample; iii) 
determining the dependent variable; iv) determination 
of the independent variables; v) descriptive 
assessment of the variables; vi) graphic AID; vii) AID 
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evaluation table, and finally; viii) graphical 
representation of the [15]. 

The stopping rule for the division of the groups 
obeys two criteria: 

• when the next division result in a group with less 
individuals that determined (in case it was determined 
that the groups should have minimal, if this study, 2.5 
million individuals), and; 

• When there is more variable among the 
predictors that discriminate between the groups as 
determined level of significance. The choice of 
predictive variables is based on a test involving the 
chi-square statistic, and the example was set a 
significance level of 5%. 

One advantage of this algorithm is that stop the 
growth of the tree before the problem of over fitting; 
another is that the result can be interpreted graphically. 
One of the disadvantages of AID exhausting algorithm 
is that require large amount of data to be able to 
ensure that the amount of observations of leaf nodes is 
significant [15,16] 

3. Results and Discussion 

This paper considers data from the Census of 2010, 
though, present deficiencies such as lack of other 
variables of interest, such as whether the disability has 
been acquired at birth or after a certain age who were 
not included in the complete questionnaire research 
that has been treated more comprehensively the issue 
of persons with disabilities. These data are useful for 
identifying the number of people with visual, hearing, 
mobility and mental, as well as its corresponding 
severity level to better assess and provide support for 
better targeting of public and private resources that 
can help improve the quality of life of these people. 

In this study, were considered following variables 
allocated on the following topics: 

(1) identification: ID number, state, region, 
household status, sex, age, race, birth registration, 
birth and nationality; 

(2) Disability: permanent disabilities to see, hear, 



move and mental; 

(3) Instruction: reading and writing, attends school 
or daycare, currently attending the course, another 
graduate course that attended and highest level of 
education; 

(4) Family: the nature of marriage, marital status 
and number of children, and finally, and, by the end; 

(5) Work: income returns home from their 
employment every day, time offset between home and 
work, occupancy condition, occupancy status, 
occupation position in main job and position in 
secondary employment. 

For this study we also created the following 
variables with their respective categorizations: on the 
topic work: income in minimum wages (mw, which at 
the time of the completion of Census 2010 was 510 
reals) categorized as 1, gain of 0 to 1 mw; 2, 1 to 3 
mw; 3, 3 to 7 mw; 4, 7 to 10 mw, and, finally; 5, 10 or 
more mw; main job (1, registered workers; 2, civil 
servants and military; 3, unregistered workers; 4, 
self-employed; 5, employers; 6, unpaid worker, and, 
by the end; 7, worker in production for own 
consumption) the topic identification: Identification 
(number questionnaire or individual) and categorized 
age (1, if 0 to 15 years; 2 if you are aged 15 to 65, and, 
finally; 3, 65 or more; race(l, white; 2, black; 3, 
yellow; 4, mulatto, and, finally; 5, indigenous); in 
family topic: number of children categorized (1 for 
free son; 2 for the number of children between 1 and 2; 
3 for many children between 3 and 5, and, finally; 4, 
number of children from six or more), and, finally; the 
topic deficiency problems: deficiencies (number of 
deficiencies that each individual possesses, and ranges 
from 0 to 4) and deficl (0 is disabled and does not 
have one if you have at least one disability). 

Tables 1 and 2 show the distribution of proportion 
in (%) for the variables total, able to read and write 
and instruction level in Table 1, gender, categorical 
age and income range categorized in Table 2 for each 
level of the variable disabilities. 

In examining Table 1 shows the variable total 76. 1% 
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Table 1 Distribution proportion in (%) of disabilities for 
the variables disabilities, total, read and write and 
instruction level. 

disabilities total read and write instruction level 



yes 



no 



1 



76.1 90.4 

17.2 82.8 
5.0 68.0 
1.6 57.3 
0.1 47.8 



9.6 
17.2 

32 
42.7 
52.2 



60.2 14.9 18.7 5.7 
64.7 13.4 16.0 5.6 
81.0 8.6 7.9 2.4 
87.7 5.9 4.8 1.5 
88.9 5.1 4.4 1.5 



Table 2 Distribution in (%) of to see for the variables 
disabilities, total, read and write, instruction level, sex, 
categorical age and income class. 



disabilities 



sex categorical age 



income class 



male female 1 



1 



4 5 



0 51.3 48.7 32.664.2 3.3 27.5 57.8 10.3 2.2 2.2 

1 45.3 54.7 10.675.1 14.3 34.7 52.1 9.2 1.9 2 

2 41.9 58.1 3.4 59.6 37.1 47.6 44.1 6.1 1.1 1.1 

3 41.7 58.3 2.3 41.0 56.7 55.4 38.6 4.4 0.8 0.8 

4 45.0 55.0 5.8 47.9 46.3 57 37.2 4 0.9 0.9 

of respondents no deficiencies considered in this study; 
17.2% one; 5% two; three 1.6%, and, finally; 0.1% 
presents four types of deficiencies studied. 

Note also that the greater the number of deficiencies 
that the person has, the lower the proportion of those 
people who can read (the proportion of 90.4% of the 
group of people who do not have disabilities followed 
in decreasing rate until at 47.8% of people with 
disabilities the four), the greater the proportion of 
people attaining maximum incomplete primary level 
(up from 60.2% the proportion of people who do not 
have disabilities following at an increasing pace until it 
reaches 88.9% of people presenting four deficiencies) 
and the lower the proportion who has full fundamental 
level or more such as instruction level, which shows 
that the greater the number of shortcomings, the greater 
the difficulty these persons to reach a better level of 
education (up from 14.9% proportion with the level of 
education among fundamental level complete and 
incomplete average, 18.7% between mid-level full and 
incomplete higher, and finally, 5.7% with college 
degrees or more for people who do not have disabilities, 
followed in descending order until; ratio of 5.1% 
between fundamental level complete and incomplete 



secondary level, middle level between 4.4% complete 
and incomplete higher, and finally; proportion 1.5% 
with college degrees or more for people with the four 
disabilities). It also shows the great disparity between 
people who does not have disabilities and those with at 
least one of the deficiencies at all levels of instruction. 

The results in Table 2 show that people having at 
least one type of deficiency is most females 
constituting 54.7% of individuals having deficiencies 
one study, 58.1% two, three 58.3%, and 55.0% four. 

Regard the different age group shows that groups of 
people who have no more than two deficiencies are 
more concentrated in the age group between 15 and 65 
years with 64.2% of the group who has no disability, 
75.1% people who have a disability, and 59.6% of 
people who have two shortcomings: the group of 
people who have three deficiencies are more 
concentrated in the age above 65 years with 56.7%, and 
the group of people the four studied with disabilities 
are more concentrated in the age group between 15 to 
65 years with 47.9% followed by the age group over 65 
years with a proportion of 46.3%, it is also possible that 
the measure increases the number disabilities greater is 
the tendency of these people to become more 
concentrated in the range between 0 and 1 with a 
minimum wage rate of 27.5% in the group of people 
that has no disability and goes increasing pace until it 
reaches 57% of people with disabilities the four 
considered in the study and decreases in other income 
brackets, and also shows the great disparity between 
people who has no disability and among those with at 
least one of the deficiencies in all income brackets. 

Tables 3 and 4 show the distributions of ratios in (%) 
for variables read and write and instruction level in 
Table 3, gender, age and income groups categorized in 
Table 4 for each level of the variable to see. 

From the results of Table 3 show that 0.2% of 
respondents can't see any way 3.3% can see with great 
difficulty, 15.1% can see with some difficulty, and 
finally 8 1 .4% do not exhibit any deficiency in visual 
terms. 
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Table 3 Distribution in (%) of to see for the variables see, 
total, read and write and instruction level. 



to see 


total 


read and write 


instruction level 








yes 


no 


1 


2 


3 


4 


0 


0.2 


69.8 


31.2 


71.8 


10.6 


12.2 


4.9 


1 


3.3 


70.3 


29.7 


77.2 


9.9 


9.9 


2.8 


2 


15.1 


81.4 


18.6 


68.5 


12.7 


15.1 


5.4 


3 


81.4 


89.2 


10.8 


62.5 


14.2 


17.5 


5.5 



Table 4 Distribution in (%) of to see for the variables see, 
sex, categorical age and income class. 



to see 




sex 


categorical aee 


income class 




male 




female 


' r r~ — 

1 2 3 


1 


f 

2 3 


4 


r 

5 


0 


47.7 


52.3 


13.6 55.8 30.6 


25.2 


56.4 12.1 


2.9 


3.4 


1 


45.3 


54.7 


5.5 63.9 30.6 


45.9 


46.0 5.9 


1.1 


1.0 


2 


41.9 


58.1 


8.0 73.4 18.6 


36.0 


51.1 9.1 


1.9 


2.0 


3 


41.7 


58.3 


31.3 64.1 4.7 


28.1 


57.4 10.2 


2.1 


2 2 



It is also observed that the higher the degree of 
severity of the disability, the lower the proportion of 
people who can read and write following 89.2% of 
people who are not visually impaired followed in 
decreasing pace to reach 69.8% in people who can't 
see at all. 

With regard to the level of instruction is noted that 
the greater the severity of visual impairment less is the 
possibility of obtaining a better educated, except that 
the group of people who can't see any way tend to 
have a better educated than the group of people who 
can see with difficulty, since the group of people who 
can not see any way has a higher proportion of people 
with education level at most elementary education and 
a lower proportion of people with complete primary 
level or more. 

Table 4 generally shows a predominance of females 
in all different severities of visual impairment. 
Regarding the age groups shows a prevalence of 
between 15 and 65 years, with people aged over 65 
years present in groups that can't see any way and that 
presents a lot of difficulty with the proportion of 
30.6% in both cases. 

Table 4 also shows that the group of people who 
can't see any way they focus more on income range 
between one and three minimum wages in proportion 
of 56.4%, the group of people who see with difficulty 
are more concentrated in the income range between 



one and three with the proportion of 46.0% followed 
by income range between zero and with a minimum 
wage proportion of 45.9%, and finally, as the group of 
people who see with some difficulty are more 
concentrated in the income range between one and 
three minimum wages in proportion of 51.1%. 

The Tables 5 and 6 present proportions in (%) 
variables total, read and write and instruction level in 
Table 5, and; gender, age and income group 
categorized in Table 6 by level to listen. 

From the results of Table 5 we can see that 0.2% of 
respondents can not hear at all, 1.0% can hear with 
great difficulty, 4.1% can hear with some difficulty 
and 94.8% did not present difficulties in terms of 
hearing. Note also that the more severe the disability 
hearing decreases the proportion of people who can 
read and write (the proportion of 88.2% of people who 
do not have hearing disabilities following in 
decreasing rate until reaching 62.9% group of people 
who can not hear at all), on the other hand, the set of 
people who can hear with difficulty find greater 
difficulty in obtaining a level of instruction rather than 
the group of people who fail in any way, by contain 
higher proportions of people with education level at 
most incomplete primary and lower proportions for 
other levels of education. 

Table 5 Distribution in (%) of to listen for the variables 
total, read and write and instruction level. 

to listen total read and write instruction level 



yes 



no 



1 



0 0.2 62.9 37.1 76.3 9.6 10.2 3.5 0.5 

1 1.0 63.8 36.2 82.6 7.8 7.4 2.2 0.1 

2 4.1 73.1 26.9 75.610.1 10.7 3.5 0.2 

3 94.8 88.2 11.8 61.714.417.9 5.6 0.5 

Table 6 Distribution in (%) of to listen for the variables sex, 
categorical age and income class. 



to listen 



sex categorical age 



income class 



male female 1 



1 



0 51.0 49.0 16.4 65.4 18.2 31.2 55.2 9.1 2.2 2.3 

1 53.7 46.3 5.5 45.8 48.7 44.8 46.4 6.5 1.1 1.2 

2 51.0 49.0 7.0 56.4 36.6 40.2 48.6 8.0 1.6 1.7 

3 49.5 50.5 27.9 66.1 6.0 29.6 56.410.0 2.1 2.1 
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In Table 6 shows that most people with hearing loss 
are males. In terms of age, groups of people who can't 
hear at all and who do not have hearing problems are 
more concentrated in the range between 15 and 65 
years, with the proportion of 65.4% and 66.1% 
respectively while groups of people who have little 
and very difficult to hear are concentrated in the 
interval between 15 and 65 years, with the proportion 
of 45.8% and 56.8% respectively, and older than 65 
years with the proportion of 48.7% and 36.6% 
respectively. 

Table 6 also shows that the group of people who 
can not hear at all can achieve better income levels 
than the group of people who can hear with great 
Difficulty, for this first group gets larger proportion of 
people with income class more than minimum wage 
and the last is more concentrated in the income class 
between 0 and 1 minimum wage with proportion of 
44.8%. 

Tables 7 and 8 show proportion in (%) of variables 
total, read and write and instruction level in Table 7, 
gender, age and categorized in income (income 
bracket) in table 8 per level to move. 

As for getting in Table 7 shows that 0.4% fail in 
any way, 2% can only with great difficulty, 4.7% can 
with some difficulty and 92.9% has no mobility 
problems. In all different degrees of severity score 
was greater predominance of educated incomplete at 
most fundamental and that the proportion of people 
who can read and write decreases as it increases the 
severity of to move disability (based on the proportion 
of 88 9% of the group of people that has no trouble 
walking, followed in descending order until you reach 
55% of people who can't move at all). Note also that 
the group of people who can hear with difficulty is 
more difficult to obtain a better level of education than 
the group of people who can't move at all, because the 
former has smaller proportions of people with level 
education teaching complete foundation or higher. 

Table 8 also shows that the problems are getting 
larger in females aged over 15 years, while the group 



that presents no problem for to move are males aged 
between 15 and 65 years. 

In Table 8, we can see that the groups of people 
who can't move in any way and do not present any 
difficulty are more concentrated in the income level 
between one and three minimum wages, while the 
groups of people who able to move with little 
difficulty and a lot are more concentrated in the 
income below the minimum wage. 

The Tables 9 and 10 present distribution of 
proportion in (%) of variables total, read and write and 
instruction level in Table 9, and; gender, age 
categorized and income group in Table 10 by level of 
mental. 

Table 7 Distribution in (%) of to move for the variables 
total, read and write and instruction level. 



to move 


total 


read and write 


instruction level 






yes 


no 


1 2 3 4 5 


0 


0.4 


55.0 


45.0 


84.1 8.3 7.0 2.5 0.2 


1 


2.0 


62.8 


37.2 


84.6 7.3 6.2 1.8 0.1 


2 


4.7 


70.0 


30.0 


79.3 9.1 8.7 2.8 0.1 


3 


92.9 


88.9 


11.1 


61.014.618.25.7 0.5 



Table 8 Distribution in (%) of to move for the variables sex, 
categorical age and income class. 



to move 




sex 


categorical age 


income class 






male 


female 


1 2 3 


12 3 4 


5 


0 


47.2 


52.8 


17.5 38.2 44.2 


26,4 55.4 12 2.8 


3.4 


1 


38.6 


61.4 


2.6 49.5 48.0 


51 42.6 4.9 0.8 


0.8 


2 


38.4 


61.6 


3.2 58.5 38.3 


47,3 44.4 6 1.1 


1.1 


3 


50.5 


49.5 


28.6 66.3 5.1 


29 56.6 10 2.1 


2.2 



Table 9 Distribution in (%) of mental for the variables 
total, read and write and instruction level. 



mental 


total 


read and write 


instruction level 






yes no 


1 2 3 4 5 


0 


1.4 


48.7 51.3 


86.6 6.6 5.4 1.3 0.2 


1 


98.6 


87.8 12.2 


62.1 14.317.65.5 0.4 



Table 10 Distribution in (%) of mental for the variables 
sex, categorical age and income class. 



mental 




sex 


categorical 


age 


income cless 




male 


female 


1 2 


3 


1 2 3 4 5 


0 
1 


54.6 
49.6 


45.4 
50.4 


16.6 68.3 
27.0 65.4 


15.0 
7.6 


54.8 39.2 4.4 0.8 0.8 
29.1 56.1 9.9 2.1 2.1 
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From Table 9 it can be noted that 1.4% of 
respondents have an intellectual disability with the 
proportion of 51.3% can't read and write and has 
educated up with incomplete basic ratio of 86.6%. 

In Table 10, shows that the majority of people with 
intellectual problems are male, with the proportion of 
54.6% and aged between 15 and 65 years, with the 
proportion of 68.3%. 

In addition, in Table 10, it appears that the majority 
of people with mental disabilities are more 
concentrated in the income level below one minimum 
wage ratio of 54.8%. 

Making a comparison between the results of a study 
of the distribution of the proportions by deficiencies in 
levels see (Tables 3 and 4), listen (Tables 5 and 6) 
around (Tables 7 and 8) and mental (Tables 9 and 10) 
it was possible verify that: 

The deficiency present in most people surveyed is 
the visual 18.6%, whereas this is less is the mental 
with only 1 .4% of the people surveyed. 

How to read and write a simple note the group with 
the worst result was formed by intellectual disability, 
whereas the group that presented the best result was 
the group formed by the visually impaired and the 
result was around 70% to people who can not see at 
all. 

Note that the groups of people who can't see, listen 
or move in any way get better educated and income 
than the group of people who can see, listen or move 
with very difficult, experts believe that this type of 
result may be due to actions of the government since 
the days of empire invested more in equipment and 
infrastructure for teaching Braille than equipment for 
low vision as telelupa, electronic magnifier, more 
wheel chair that physiotherapy treatment for can 
acquire more mobility and listening aid than treatment 
with speech therapy. 

In order to evaluate the different degrees of severity 
yeah, can not in any way but with much difficulty, 
homogeneity tests were performed for these severity 
levels of the variables see, move and hear that in all 



cases were obtained level significance of 0.000 for 
this test and concluded that there is no homogeneity 
between these levels tested. 

With this results checking that people who can't in 
any way has better performance in terms of level of 
education and income than people who answered yes, 
but with much difficulty and experts consulted this 
type of situation occurs due to public actions as 
establishing special schools, production and import of 
equipment as regulate, typewriters and braille printers, 
these products, which benefit else who can't in any 
way that anyone can, but with great difficulty, as in 
the case of the latter, often need special equipment and 
appliances just existing in the Brazilian market as 
telelupa low vision in the case of vision, devices that 
aid in hearing in case of trouble hearing and 
wheelchairs for people who have difficulty walking 
reflecting a serious economic problem for families and 
for the country itself, since these equipment are 
imported and besides the high price also adds import 
fees which are also high, and in most cases, are not 
manufactured in Brazil, and notices the nonexistence 
of public actions by various governments that best 
makes possible this kind of production and 
procurement. 



The Tables 1 la and 1 lb show distribution in 



of 



proportions of the formal labor represented by the sum 
among proportions of the registered workers, civil 
servants and military and employers, and, informal 
work represented by the sum of proportions among 
unregistered workers, self-employed, unpaid worker 
and worker in production for own consumption. 

From Tables 11a and 1 lb are possible to verify the 
predominance of informal work on the formal in all 
situations, being more pronounced in groups with 
intellectual disability (71.3%) who see (67.1%), 
listening (67.4 %) or moving (71.6%) with difficulty or 
having two (70.0%), three (76.0%) or four (74.4%) 
disabilities. 

The graphs of Figs. 1 to 5 show the diagrams for 
AID number of deficiencies (deficiencies) and 
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Table 11a. Distribution proportion in (%) of the formal 
labor and informal work for disabilities and to see 





formal 


informal 




formal 


informal 


disabilities 


labor 


work 


to see 


labor 


work 


0 


49.7 


50.3 


0 


51.1 


48.9 


1 


41.5 


58.5 


1 


32.9 


67.1 


2 


30.0 


70.0 


2 


41.1 


58.9 


3 


24.0 


76.0 


3 


49.2 


50.8 


4 


25.6 


74.4 









Table lib. Distribution proportion in (%) of the formal 
labor and informal work for to listen, to move and mental. 

formal informal fonnal informal formal informal 

to listen labor work to move labor work mental labor work 



0 48.9 51.1 

1 32.6 67.4 

2 36.0 64.0 

3 47.6 52.4 



47.3 
28.4 
30.6 
48.1 



52.7 
71.6 
69.4 
51.9 



28.7 
47.3 



71.3 
52.7 



deficiencies to see, move, mental and hearing 
respectively. For each of these diagrams were formed 
six groups: the first group (black), second (red), the 
third (blue), the fourth (green), the fifth (brown) and 
sixth (pink). 

Studying graphic diagram of Fig. 1 AID shows that 
the variables that most contribute to the appearance of 
a smaller number of disabilities are unmarried, 
economically active, attending a course and do not 
have secondary employment. 

Among the different profiles of the groups formed 
what best contributes to fewer disabilities is the group 
formed by four unmarried and attending a course, 
while the profile that contributes most to the increase 
in the number of disabilities is the group 3 formed by 
being married, divorced, widowed, legally separated 
condition and not economically active. 

For the case of Fig. 2 shows that the variables that 
most contribute to a lower incidence of visual 
impairment are unmarried, activity status, 
economically active status, level of education from 
primary level onwards and attended or attends 
daycare. 

Considering these observations, the better is the 
group formed by 6 which presents unmarried, have 
attended daycare, education level completed elementary 
onwards followed by group 4 unmarried and attends 
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public or private kindergarten, while the largest 
contributor to the increase visual impairments is group 
5 consists of being married, divorced, widowed, 
separated or legally separated and not economically 
active. 

As for the Fig. 3, the variables that most 
contribute to a lower incidence of deficiencies for 
around are unmarried, economically active condition, 
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already attended or attends public or private day care 
and returns home from work every day. 

Among the different groups formed what best 
contributes to a lower incidence of disability for the 
group was around 4 formed by unmarried and attends 



public or private daycare and that contributes to the 
worsening of this deficiency was the group formed by 
three married state , divorced, widowed and separated 
or legally separated situation and not economically 
active. 

In the case of Fig. 4, the variables that most 
contribute to a lower incidence of mental disabilities 
are economically active status, level of education from 
the complete primary onwards and female and the 
group formed that contributes to the higher incidence 
of mental problems was 4 formed by economically 
active status and education level of at most incomplete 
primary. Among the groups, that best contributes to 
the four, lower incidence of mental or intellectual is 
the group formed by six economically active status, 
level of education completed elementary or older and 
female. 

Finally, in the case of Fig. 5, the variables that most 
contribute to a lower incidence of hearing impairment 
are unmarried, economically active condition, attends 
course from regular elementary and fundamental 
education level or more. 

The group that best contributes to a lower incidence 
of hearing impairment is the group fonned by four 
unmarried and attends regular course of school or 
more and profile that contributes to the increased 
incidence of hearing impairment is the group formed 3 
by marital status married, widowed, divorced, 
separated or legally separated and not economically 
active situation followed by group 5 formed by 
unmarried, attends literacy course for young adults 
and educated incomplete primary or no education. 

The variables that were selected for this study were 
provided activity in the five cases studied: marital 
status in four; education level in three; course and 
attends daycare in two, and finally, secondary job, 
returning home from work, sex and nature of the 
union are present in every single situation. 

By doing an analysis of the plots in Figures 1-5 we 
can see that except mental disability (Figure 4), the 
first variable that partitions the other dependent 
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variables are different marital status, which in turn 
discriminates economically active status variable. 

4. Conclusions 

The greater the number of disabilities tends to be the 
lowest level of education and income that these people 
will get in their work. 

The level yes, I can't in any way tends to achieve 
better levels of education and income level rather than 
the level, but very difficult for variables to see, move 
and listen. 

Note the existence of a large gap between people 
with disabilities and without disabilities as the level of 
education and income level, with the highest 
concentration of people with disabilities in the income 
ranges between 0 and 1 mw and between 1 and 3 mw, 
and; educated to uneducated incomplete primary. The 
higher the level of education and a higher income 
bracket, the greater the inequality between disabled and 
non-disabled. 

The profile that most contributes to the incidence of 
these deficiencies is formed by being married, 
widowed, divorced, separated situation legally and 
economically inactive. 

A deficiency in this population is the most visual and 
lower this is mental. 

Deficiency that has the highest amount of people 
who can't read and write is the mind while the one with 
the greatest amount of people who can read and write is 
visual. 

Group of people who see, hear or move achieve 
results with some difficulty, although smaller, but very 
close to people who do not show any deficiency in 
terms of level of education and income. 

The value of this research lies in the fact that there is 
virtually no work to do comparative study between 
people with disabilities and who do not have 
disabilities considering explanatory variables such as 
education, income, age and others. 

As a continuation of this work in the future we 
propose: 



First: new study comparing people who belong to 
groups that can not see, hear or move in any way they 
see and with some difficulty; 

Second: Evaluate together with the possibility of 
improvements IBGE's census questionnaire adding the 
topic disability for each different disabilities a 
question that asks whether the disability is from birth 
or after birth was that age, justified by the fact people 
who acquired disability after a certain time the most 
difficulties to adapt this new condition which carries a 
certain time to adapt in their new condition; 

Third: Application of new analysis using other 
techniques such as geostatistics, canonical correlation, 
structural equation modeling and other that may assist 
in the investigation between different disabilities and 
other variables considered as predictors; 

Fourth: periodic repetition of the survey data in 
order to allow better monitoring of people with 
disabilities and allow other comparative studies can be 
made; 

Fifth: Consider the relationship between the 
different disabilities and other variables collected in 
the census as a race, and other major work; 

Sixth: Consider also the study by region, state and 
county, and finally; 

Seventh: Submit new studies and monitoring of 
variables from the categories that contribute in 
enhancing the various deficiencies after using the 
technique AID as analysis. 

The results of this work can benefit public 
managers in better support in the care of persons with 
disabilities to have knowledge of who they are, where 
they are and how. 
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Abstract: Banking institutions all over the world face significant challenge due to the cumulative loss due to defaults of borrowers of 
different types of loans. The cumulative default loss built up over a period of time could wipe out the capital cushion of the banks. The 
aim of this paper is to help the banks to forecast the cumulative loss and its volatility. Defaulting amounts are random and defaults occur 
at random instants of time. A non Markovian time dependent random point process is used to model the cumulative loss. The expected 
loss and volatility are evaluated analytically. They are functions of probability of default, probability of loss amount, recovery rate and 
time. Probability of default being the important contributor is evaluated using Hidden Markov modeling. Numerical results obtained 
validate the model. 

Keywords: Random point process, expected cumulative loss, non Markovian, hidden Markov model. 



1. Introduction 

For most of the banks, loans are the largest and the 
most obvious source of credit risk. The banks should 
ensure adequate capital available against this risk. The 
higher the risk more chances of default and resultant 
losses to bank. Hence it is very important for banks to 
measure, monitor and control credit risk over a period 
of time. Banks all over the world face significant 
challenge due to default of borrowers. The turbulence 
that dominated financial markets throughout the world 
during the last three decades has proved the 
ineffectiveness of the standard risk management tools 
used. Most of these tools are based on the normal 
distribution and they do not account for the random 
events. Therefore a better modeling tool should be 
based on random point processes because defaults 
occur at random instants of time and defaulting 
amounts are random. 

The cumulative loss due to defaults occurring at 
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random instants of time and built up over a period of 
time could affect the income of banks and could reduce 
the capital cushion. Though large number of research 
publications dealing with default emphasize on 
Markovian approach, authors like Nickell, Peraudin 
and Varetto (2000), Bangia et al (2002) and Coudere F 
(2005) insist on including non Markovian features of 
the problem of cumulative loss. Agencies like Moodys 
and Standard & Poor have refined their initial ratings 
by adding certain variants to accommodate the non 
Markovian features of the problem. In this paper a time 
dependent non Markovian random point process is 
studied to model the expected cumulative loss and its 
volatility. 

The aim of the study is to derive an compact 
expression for expected cumulative loss and volatility 
so that banks can use them with the data available. 
Random point process form a special class of stochastic 
process used in modeling wide range of problems with 
the following common characteristics. 

Events occur at random instants of time and each 
event triggers a feature like random default loss. The 
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quantity that is of interest is the cumulative nature of 
the feature over a period of time. Stochastic processes 
associated with random points on a line were studied 
by Ramakrishnan (1953) and subsequently studied by 
S.K.Srinivasan and KSS Iyer (1970). In the present 
paper it is proposed to consider that each default occurs 
at random instant of time, not necessarily of Poisson 
type, triggers a random loss and the amount builds up 
over a period of time. The expected cumulative loss 
and its volatility are explicitly obtained and they are 
functions of probability of default, expected loss 
amount, recovery rate and time. 

The paper is organized as follows: In section 2 
random point process model is presented for the 
cumulative default loss. Section 3 considers the 
formulation of stochastic differential equation satisfied 
by the probability frequency function of cumulative 
loss. Using Fourier transform closed form analytic 
solutions are obtained for the first two moments of the 
cumulative loss for the Poisson arrival of defaults. For 
non Poissonian arrival a hidden Markov modeling is 
used to derive the probability of default in section 4. 
The maximum likelihood probability of default among 
borrowers of the same type of loan by classifying them 
into three states like least risk, medium risk and high 
risk. In the absence of real data the numerical results 
are obtained in terms of dimensionless ratio and the 
model is validated by the graphs. 

2. Random Point Process Model 

Consider defaults of type j occurring at random instants 
of time ti, t 2 , .... t n . We associate with each tj a random 
loss a J (t ; ), where a J (t ; ) are statistically independent 
random variables. We associate a time dependent 
deterministic hazard function h J (t, tj) to capture the 
impact of the unit default loss of type j at t ; at a later 
time t where tj< t. This could represent the built up rate 
of loss function over a period of time. For a typical 
realization of defaults at ti for i = 1, 2, 3. . . N, the 
cumulative loss at time t may be written as 

L J (t) = £ a J (ti). h J (t, ti) for i = 1, 2, (1) 



a j (t 2 ) 



a j (ti) 



a j (t 3 ) 



U(t) = lai(t i ).hi.(t-t ] ) 



tl 

Fig.l 



t2 



t3 



It may be noted L J (t) is a random variable because 
a J (t ; ) is random. Also L J (t) is non markovian and the 
normal markovian analysis cannot be carried out. 
However the probability frequency function for L J (t) 
can readily give the moments of L J (t). 

2L j (t) for j = 1, 2, 3....k will give the total 
cumulative loss due to K types of loans. For the rest of 
the discussion, we drop the index j, as the solution for 
L(t) will be identical except for the default and 
recovery rates . The first moment of L(t) gives the 
expected cumulative loss at time t and second moment 
gives the volatility at time t. 

3. Derivation of Stochastic Differential 
Equaton for L(t) for Poisson Arrival Rate of 
Defaults. 

In this section we obtain the stochastic differential 
equation satisfied by Ji(L,t) the probability frequency 
function governing the cumulative loss function L(t). 
We assume for the present discussion arrival of 
defaults obey Poisson distribution though the results 
obtained are equally valid for non Poissonian 
distribution for deterministic time dependent default 
rates. 7i(L,t)dL denotes the probability that L takes a 
value between L and L+dL at time t assuming L to be 0 
at time t = 0. With each default randomly occurring at tj, 
we associate a loss a(t ; ). h.(t - tj) where probability of 
random loss a(tj) is given by p(a ; ). 

Analyzing the possible events in the time interval (0, 
At) there are two possible events 

(1) A default occurs in (0, At) with probability A,. At 
in which case a contribution of a.h(t) to L(t) arises and 
the corresponding contribution to n(L, t) is given by the 
probability 
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XAt . j a 7t(L - a.h(t), t). p(a) e^' da 

where the exponential term ensures that the probability 
of not coming out of default during the time t and |i 
represents the recovery rate. p(a) represents the 
probability of loss amount is a. 

(2) No default occurs in (0, At) with probability (1 - 
XAt) in which case the corresponding contribution to 
7i(L, t) is 

(1-AAt). n(L, t-At) 

Combining the above two we have partial 
differential equation for n(L, t) 

dn(L,t)/ St = -Xn(L,t) + Xj p(a). 

n(L - ah(t), t ) e" 1 ". da (3.1) 

The above equation has been solved by using Fourier 
transform and by using the properties of characteristic 
function the moments of L(t) are derived. 

Solution 

l 00 

Defining n (s,t) = — J 7i(L,t) e lsL dL equation 3.1 

o 

takes the form 

dx (s,t)/a = 

00 

-X„ (s,t) + \x^ (s,t).p(a)e iash(t) e^ da (3 .2) 
o 

Equation (3.2) gives the characteristic function of 
L(t) viz 7i(s,t). 

This yields all the moments of L(t). 

Defining (p(u) = Jp(a) e lua da and substituting in (3.2) 
we get 

3 x (s,t)/dt = —%x (s,t) + Xn (s,t)e-^(p(sh) (3.3) 
The solution of the above equation can be written as 

t 

n (s,t) = exp [ -Xt + X J e" 1 " <£(sh)dx ] (3.4) 

o 

By taking inverse Fourier transform of the above 
equation we can obtain n(L,t). 

00 ( 

jc(L,t) = ^ J exp[-isL + X^f J A.h"(T)e" llI E[a"(T)]dT]ds (3.5) 
0 0 



The n moment of the cumulative loss L(t) is given 

by 

E[ L(t)] n = i n .5 n /5s n . {7r(s,t)}| s=0 (3.6) 
Equation (3.5) gives the probability of cumulative 
loss L at time t. 

The expected cumulative loss and the volatility can 
be obtained from the first two moments of L(t). 
Differentiating Equation (3.3) gives the moments of L(t) 
and the first two moments are the first and second 
derivatives of n (s,t) with respect to s and for s = 0. 
Denoting N(t) and M(t) for the first two moments of L(t) 
we get the following equations satisfied by N(t) and 
M(t). 

aN/at=-^N+^Ne _,it O(0)+h(t)O'(0) Xe^ (3.7) 
dMldt= 

-XM+Xe'^ [M<£(0)+2Nh(t) 0>'(0)+h 2 (t) 0>"(0)] (3.8) 
Solving (3.7) equation we get, 

N(t)=h <D'(0) Xe m Je^^dT (3.9) 

o 

Where 

> P(t) = lt-^(0)/|i(l-e" llt ) 
O(0)=p 0 (a m -a 0 ) 
0)'(0)=po (a m 2 -a 0 2 )/2 
a m =Maximum Loss amount 
a 0 =Minimum Loss amount 
Solving (3.8) equation we get, 

M(t) = J ' e x " 12/2 0(0). Xe ,Jt [2N(t) h 2 (t) 3>'(0)+ h 2 (t) 

o 

O"(0)]dt (3.10) 
M(t) gives the second moment and the volatility of 
the cumulative loss can be evaluated from M(t) and 
N(t). It may be noted it depends on X, |i, a m , a 0 ,po and 
h(t) the built up rate of loan over a period of time t. 

The expected cumulative loss depends on 
probability of default X, recovery rate (X, built up rate of 
defaulted amount h and po the probability of the loss 
amount. The probability of default is an important 
measure for banks. In the next section, we discuss the 
hidden markov modeling approach to evaluate the 
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probability of default. 

4. Hidden Markov Modeling to Evaluate 
Probability of Default 

Hidden Markov models have been used in variety of 
applications like Speech Recognition, signal processing, 
artificial intelligence, computational biology, finance, 
image processing and bio statistics. These models 
possess the ability to present unobservable 
developments that take place with observable out puts. 
A comprehensive treatment of Hidden Markov models 
can be found in Rabiner (1989). 

In this section we propose to evaluate probability of 
default using hidden markov model. Hidden Markov 
modeling consists of two sets of states namely hidden 
states and observable states. Hidden states are 
described by Markov Process and their transition 
probabilities are given by transition matrix. 

The observed output stochastically depends on 
hidden states. Emission probabilities from hidden 
states to observable states are given by the emission 
matrix. The initial probabilities for the hidden states are 
known. 

Elements of a discrete HMM are: 

(1) Hidden states Q = [ q ; ], i = 1, 2,. . .N 

(2) Transition probabilities A = { ay = P(qj)at t+1/ 
P(qO at t}, where P(a/b) is the conditional probability of 
a given b and t = 1, 2, ... T. 

(3) Observations O = {o k }, k= 1,2,... M 

(4) Emission probabilities B = { bjk = bj(Ok) = 
P(Ok/qO } . B is the probability that the output is Ok given 
the current state is q;. 

(5) Initial state probabilities n = { p; = P( q; at t = 1) } 
The model is characterized by the complete set of 

parameters X = { A, B, n } 

The three important problems that can be solved 
using HMM are: (a) given the model the computation 
of the probability of the occurrence of the observation 
sequence (b) given the model the sequence of the states 
that maximizes the joint probability of state sequences 
and observations sequence (c) adjustment of the 



parameters of the model to maximize either of the 
above two probabilities. 

5. Working Methodology for Calculating 
Probability of Default 

We propose to consider the population of borrowers 
of type j loan in three groups like least risk, medium 
risk and high risk based on the risk levels from past 
data of the borrowers available with banks. We 
consider the number of states as four including the 
default state. State 1 is least risky, state 2 is medium 
risk, state 3 is high risk and state 4 is default state or the 
absorbing state. The first three states generate the 
observable states namely risk levels and the state 
transition takes place from left to right among the states 
with default as absorbing state as shown in Fig. 2. 

We assume that default can come from any of the 
three groups of borrowers with three different risk 
levels. Default is the absorbing state and the best 
sequence is the one that maximizes the probability for 
the given observation sequence ending with default 
state. For simplicity it is assumed that the transitions 
among hidden states are considered from left to right. 
In other words the transition matrix is not Ergodic. 
Viterbi's algorithm chooses the best state sequence that 
maximizes the likely hood of the state sequence for 
given observation sequence when the Markov process 
is not ergodic. For ergodic matrix posterior decoding 
can be used. 

Fig. 3 depicts numerical results from Viterbi 
decoding of the hidden markov model of Figure 2. The 
plot shows variation in total probability of default with 
increasing risk levels of the low, medium and high risk 
states. The levels of risk are implemented as the 
emission probabilities of the observable states Oi-C>6. 
Two sets of risk levels have been studied to ascertain 
the variation in default probability with risk. From the 
plot it is clear that the probability of default increases 
gradually with increasing risk levels of the various 
states in the model. 
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Numerical Results: In the absence of real data, to 
gain some insight into the validity of the model, we 
have plotted the expected cumulative loss in terms of a 
non-dimensional ratio of maximum loss to minimum 
loss for default. N(t) is the expected cumulative loss 
and is function of X,[i,h(t) and p(a). Assuming h(t) and 
p(a) as constants h and p 0 we get , 

N(t) = ^p 0 h(r 2 -l)/2 . ( (aoV- e^V]) / ([ 
3>(0))-fi])) 

Where r = a m /a 0 

The above expression can be written to the first 
approximation in a non-dimensional form 

N(t)/(^p 0 hao 2 ) = ((r 2 -l)/2).t 

Fig. 4 represents the above graph for different values 
of t and r. 



The graph indicates at t=0 the cumulative loss is zero 
for all values of r. Also when r =1 namely a m = a 0 the 
cumulative loss is zero. The cumulative loss increases 
with increasing time and values of r. 

Time required for cumulative loss to reach 
critical value. 

Assuming the cumulative loss N(t) reaches the 
maximum value N c after a time T c . we have plotted T c 
for different values of from the following expression 

T c = ((2N c )/(r 2 -l)).(l/(^p 0 ha 0 2 .l)) 

From the graph in Fig. 5 it is seen that the time taken 
to reach critical cumulative loss falls rapidly with 
increasing values of r and the curve is asymptotic for r 
= 1. 
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6. Conclusion 

In this paper an attempt is made to model the 
cumulative default loss for different types of loans as a 
random point process with defaults occurring at 
random instants of time. The expected cumulative loss 
for a homogeneous type of loan at any time is 
expressed in terms of probability of default, recovery 
rate, built up rate of loss amount over a period of time 
and the probability of loss amount. The total sum of the 
cumulative loss due to all types of loan can be 
evaluated from the expression £E{L J (t)} for different 



values of j = 1, 2, ... etc. We presume banks can 
calculate the above expression from the data available 
to them. The time dependent cumulative default loss is 
plotted against the non-dimensional ratio of maximum 
loss to minimum loss for default in the absence of real 
data. The time to reach the critical default loss is also 
plotted. The expression for the second moment of the 
cumulative loss, which will give the volatility, is also 
evaluated. Results obtained in this paper are equally 
valid in more general case when the probability of 
default depends on time. It is expected that the analytic 
solution derived in this paper could provide a rational 
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basis for banks to forecast the time to cumulative 
default loss that could wipe out the capital. 
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Abstract: Data are critically important to good decision-making. In an increasingly complex economy, conventional data collecting 
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1. Introduction 

Economies are constantly faced with new challenges. 
To remain relevant, official statistics have to keep up 
with the rapid changes of modern times, which 
typically require the availability of commensurate 
statistical data that users may exploit in an accurate and 
reliable way. Policy-makers, financial supervisors and 
regulators, just to name a few, require as much rich and 
timely information as possible to take appropriate 
decisions. 

The Banco de Portugal (hereinafter referred as "the 
Bank") - or any other major producer of official 
statistics, for that matter - has to ensure that the 
statistics for which the Bank is accountable retain 
relevancy over time and are able to cope with the speed 
and the scope of the main stakeholders' ever-increasing 
demand for comprehensive, detailed and high-quality 
information. 

However, the process of continuously adapting the 
statistical output to new phenomena has a number of 
serious limitations. Conventional data collecting 
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systems cannot simply keep on expanding indefinitely 
to cope with the ever-increasing need to fill the 
information gaps perceived by the users or in 
anticipation to their possible future data requirements. 
Amongst the possible motives for not pursuing 
recurrently this approach one could point out, inter alia, 
the following: 

• The resulting overburdening of respondents goes 
against well-established best practices. 

• The related initial and maintenance costs are far 
from being negligible, both to the agency that collects 
the data and to the respondents. 

• New statistical datasets (or significant 
enhancements to existing ones) require lengthy 
preparation time (years, rather than months) and, once 
launched, are supposed to remain in operation for a 
prolonged period of time (typically around five years, 
in the case of the Eurosystem statistical reporting 
systems). This time-lag could even be further extended, 
should the revision result from a major methodological 
change, as it is often the case. 

• Ad hoc surveys are, in general, too 
time-consuming and expensive, not to mention reliant 
upon the willingness to participate on the part of the 
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target population. 

In fact, the response given by conventional data 
collecting systems to new statistical demands - 
stemming from, e.g., the need to conduct 
macro-prudential analyses or to accommodate new 
data requirements related to the Bank's participation in 
the European System of Central Banks (ESCB) - is 
problematic, costly and could possibly turn out to be 
counterproductive, which helps to understand why 
more and more central banks have been reusing the 
available micro-data, thus recognising that such 
information is both useful and necessary to respond to 
the data requirements of the complex world we live in, 
and to better address new issues and challenges as they 
arise. 

Throughout this paper the term "micro-data" will be 
used to refer data about individual persons, households, 
businesses or other entities; it may be data directly 
collected by the Bank or obtained from other sources, 
such as administrative sources. 

2. From Direct Integration to Integration via 
a Data Warehouse 

2.1 Changing the Statistical Paradigm 

Data are critically important in making 
well-informed decisions. Poor quality data or, a fortiori, 
lack of data can lead to an inefficient allocation of 
resources and imposes high costs on the society. In 
ever more complex economies, the traditional 
approach to the compilation of official statistics - i.e., 
producing standard statistical tables that can only 
address a set of predefined questions - is becoming 
increasingly insufficient and ineffective. 

The Bank's strategy to deal with the challenge of 
maintaining its statistics relevant to the users in a 
shifting and more demanding environment, while 
attending to the need to keep the reporting burden on 
respondents at an acceptable level, was to enhance the 
overall efficiency of the statistical framework by 
further exploring the largely unused statistical potential 
of already existing data sources. In fact, statistically 



edited micro-data, which include e.g. data from 
administrative sources not originally intended for 
statistical purposes or even data related to the Bank's 
prudential supervision function, offer an unusual array 
of interesting features, inter alia: 

• Very good coverage of the population in most of 
the cases. 

• Relatively low reporting costs, thus helping to 
mitigate the constraints imposed by the response 
burden of the reporting agents. 

• Increased flexibility and agility as regards the 
compilation of new statistics, e.g. related to financial 
and other structural innovations. 

• More rapid response to ad hoc data requirements 
from the users - in many cases, almost in real time. 

Moreover, the evolution in network and 
communication protocols, database systems and 
multidimensional analytical systems has somewhat 
removed the potential disadvantages of having to deal 
with the huge amounts of data normally associated with 
the handling of micro-databases. (Aguiar et al, 2011) 

Best practices in compiling official statistics 
advocate that all data should be collected only once: 
any form of double reporting or redundant collection 
should be avoided and, if existing, be terminated. 
Accordingly, data already available - due to whatever 
reasons - should be reused, if found useful, for 
statistical purposes. Obvious candidates are data from 
existing Central Credit Registers, as well as data from 
Central Balance Sheet Offices databases and 
information collected within the framework of the 
Bank's prudential supervision function. The 
experience of the Bank in this area has shown that the 
use of such information for statistical purposes can lead 
to a significant reduction of the response burden, 
higher data quality and lower costs. 

On the national level, a formal exchange of 
administrative data with institutions outside the central 
bank, like the national statistical institute (NSI) or the 
tax authorities, would also help to reduce the reporting 
costs. An important precondition would be the 
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maintenance of common company registers with the 
NSI. Extending this idea across national borders, one 
could think of common international databases - e.g., 
exchanging micro-data on significant cross-border 
mergers and acquisitions that need to be recorded 
symmetrically in the respective statistics of both 
affected countries. (Liebscher et al, 2008) 

2.2 Micro-Databases Managed by the Banco 
de Portugal 

For the last 10 years the Bank has been developing 
and maintaining several micro-databases based on 
item-by-item reporting and has been exploring the 
statistical potential of these complementary sources of 
information with significant positive impacts on the 
overall quality of its statistical output. 

The databases managed by the Bank's Statistics 
Department include: 

• The Securities Statistics Integrated System (SSIS) 
database, a security-by-security and an 
investor-by-investor database that provides, in a single 
repository, data on the securities issues and holdings 
required by the different statistical domains (e.g., 
monetary and financial statistics, external statistics, 
securities statistics and financial accounts), thus 
replacing the separate and distinctive data storing 
systems that were previously in place. 

• The Central Credit Register (CCR), an 
administrative database that stores credit-related 
information supplied by all the resident credit-granting 
financial institutions. 

• The Central Balance Sheet Database (CBSD), 
which stores granular information on virtually all the 
resident corporations, collected through the so-called 
Informagao Empresarial Simplificada (IES), a joint 
effort of four distinct Portuguese public entities - the 
Ministry of Finance, the Ministry of Justice, Instituto 
Nacional de Estatistica (the Portuguese NSI) and the 
Banco de Portugal - consisting of yearly submissions 
of information by corporations, in a single, paper-free, 
electronic form, to fulfil reporting obligations of 



accounting, fiscal and statistical nature. 

Besides complementing and helping to cross-check 
the information gathered through the conventional 
channels, these micro-data have proved to be of great 
importance to the understanding of the developments 
in the Portuguese financial system, especially in the 
wake of the recent global financial crisis. 

So far, this approach has permitted, inter alia: 

• Improving the responsiveness to new users' 
requirements, particularly those arising from ad hoc 
information requests, with proven results in reducing or 
eliminating data gaps and in monitoring and assessing 
the evolution of the Portuguese financial system. 

• Curtailing the follow-up procedures as regards 
data collecting schemes, whereby respondents are 
re-contacted after the initial submission of data, to 
obtain missing information and/or to verify and, if 
necessary, to correct questionable data. 

• Enhancing the quality control procedures (e.g., 
by cross-checking elementary/raw data from different 
statistical domains), thus increasing the efficiency of 
the production process and improving the quality of 
end-products. 

• Avoiding data redundancy, while at the same time 
expanding significantly the range of statistics 
available. 

As an example, the use of the available 
micro-databases for the compilation of the Portuguese 
flow-of-funds within the national financial accounts 
has been extremely helpful, as it allows for a much 
better understanding of the interlinks within the 
resident economy and vis-a-vis the rest-of -the -world. 

2.3 Deepening Data Integration 

In keeping with such course of action, the Bank has 
been developing an approach that, once completed, 
will allow for a higher level of integration of the 
available administrative and survey data. The goal is to 
achieve a significant enhancement of the basic data 
infrastructure without jeopardizing the provisions in 
the legislation, codes of practice and protocols that 



254 



Impact and Benefits of Micro-Databases' Integration on the Statistics of the Banco de Portugal 



protect data confidentiality. In addition, a reduction in 
respondent burden and an increase in the breadth and 
depth of the information available to policy-makers 
and researchers are expected. 

An architecture based on business intelligence - 
broadly defined as a category of applications (e.g., 
decision support systems, query and reporting, online 
analytical processing, statistical analysis, forecasting, 
and data mining) and technologies for gathering, 
storing, analyzing, and providing access to data, to help 
a variety of users to make better business decisions 
(Terzic, 2008) - could significantly contribute to 
meeting the Bank's concerns in this area. With this in 
mind, the Bank set off a study in 2008 aiming at 
defining a business intelligence framework to be used 
as a reference in all future information technology 
developments in the statistical realm. This framework 
will be built upon three pillars (Aguiar et al., 2011): 

• A common technological infrastructure across the 
various information systems, to facilitate the 
integration and re -usage of components and to promote 
data access efficiency and transparency to final users. 

• A centralised reference database, to provide 
common reference data (e.g., identification criteria for 
the relevant entities that are observed, characterization 
of variables and classifications) and to enable the 
linkage of information across different sources and 
systems. 

• A data warehouse approach, to guarantee a central 
access point to all statistical data, independently of the 
input source or the production process. This implies, 
inter alia, a data structure specified on the basis of 
common criteria, valid across the different sets of data. 

At the moment, the Bank's statistical information 
subsystems are in the process of being reformulated 
according to this model: the SSIS and the balance of 
payment and international investment position 
statistics, on one hand; the CBSD and the CCR, on the 
other hand. 

Data integration is concerned with integrating unit 
record data from different administrative and/or survey 




Fig.l From direct integration to integration via a data 
warehouse. 

sources to compile new official statistics which can 
then be released in their own right. Integration of 
micro-data is a powerful approach to enrich the already 
available information - e.g., by allowing efficient 
cross-data comparisons and quality checks among the 
different statistical domains. Surely, this is easier said 
than done - it is a rather complex process and, pending 
on the degree of integration to be achieved, it can be 
characterized by different features. 

The prevailing, stand-alone, islands of information 
may be very diverse, making it technically difficult to 
create homogeneous information systems. In addition, 
there are many different practical and methodological 
problems that must be previously addressed when two 
or more sources are to be integrated, inter alia (Di Zio, 
1998): 

• Harmonising populations - e.g. determining the 
group of entities that belong to a given institutional 
sector (financial and non-financial corporations, 
general government, households and non-profit 
institutions serving households) -, identification 
criteria, reference periods (annually, quarterly,...), 
variables and classifications. 
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• Adjusting for measurement errors (accuracy) and 
for missing data. 

• Deriving variables. 

However, such shortcomings may very well be 
offset by the possible benefits of integrated data sets. 
The latter include, according to UNECE (2009): 

• Compiling new or enhanced statistics. 

• Producing more disaggregated information for 
measures where some information currently exists. 

• Carrying out research using composite micro-data 
that cover a wider range of variables for a larger 
number of units than available from any single data 
source. 

• Potentially improving or validating existing data 
sources. 

• Possibly reducing respondent burden. 

These benefits could be illustrated by the following 
case, extracted from the Bank's own statistics: a given 
corporation, providing annual accounting data under its 
IES reporting obligations, might also be answering to 
the Bank's IS II survey (Inquerito sobre Investimento 
International) and, at the same time, having its 
securities issues and holdings recorded in the SSIS 
database; in an integrated system, it would be possible 
to ensure the compatibility of these data at a 
micro-level, thus providing a powerful tool for the 
compilation of financial accounts (which require that 
total uses equal total resources in the domestic 
economy). Nonetheless, as referred above, having a 
partial integration, e.g. one that allows for a unified 
view on two different sub-systems like the CBSD and 
the CCR, clearly enriches analytical data awareness. In 
fact, the idea of pooling together all the data on 
financial or non-financial corporations available at the 
Banco de Portugal is rather appealing; it would allow 
us to have a more specific and detailed view on this 
particular institutional sector. In the case of the 
financial corporations, the advantages would be even 
stronger should we take into account the information of 
yet another important subsystem: the data used for 
supervisory purposes. 



4. Concluding Remarks 

The implementation of an architecture framework 
such as the one summarized above will contribute to 
the construction of a coherent and integrated statistical 
system as opposed to having multiple systems that 
coexist but are not connected in an efficient way. 

Such approach has only been possible because of the 
possibilities brought in by the information technology 
(IT) revolution. But even though IT has enabled the 
statistical community to carry out the current 
procedures for collecting, compiling and disseminating 
statistics more efficiently, albeit at a non-negligible 
cost, it is important to reflect on how such revolution 
can be used to introduce new and more effective 
procedures. 

Benefits are evident but there are also problems, 
challenges and cautions with the use of integrated 
micro-data, particularly those related to confidentiality 
issues. As said before, data already available should be 
reused if found useful for (other) statistical purposes; 
that being the case, it is necessary to strictly safeguard 
their confidentiality and to ensure that the sharing is 
legally allowed or explicitly agreed by the reporting 
agents. However, because of legal constraints, 
confidentiality makes the access to some useful data 
sources problematic and disclosure is a constant 
problem when we need to release data. 

A data integration process is complex and can be 
characterized by different steps. One of these steps is 
adopting a unified view on the existing micro-data data 
sources creating a customized view on a sub-set of data 
{e.g. the financial or non-financial sectors). 

Integrated micro-data have the potential to support, 
if need be, the drilling down of the most summarized 
levels of data to the most detailed ones, which may help 
to confirm (or to disprove) trends and developments 
conveyed by macroeconomic statistics and, 
concomitantly, to explore and/or to elucidate their 
possible implications for e.g. financial stability 
analysis and systemic risk assessment. (D'Aguiar etal, 
2011). 
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1. Introduction 

The object to identify features, to compare them with 
the data from the library and infer the presence 
probability of anomalies, to preallocate to a plurality of 
present on a specific image. In most cases, the 
investigated images are present noise, the texture 
similar distortion to those regions owned by the object 
under investigation. All this complicates process of 
selecting objects and correctly display their borders, so 
contouring and segmentation algorithms have a very 
important role in the automated processing. 

The operator Sobel - this is one the best algorithms 
available border selection, It is often applied as one of 
steps a more difficult and exact algorithms (such as an 
operator Kenny). 

2. Experimental Section 

The Sobel operator used in image processing. Often 
it is used in the algorithms border selection. This is a 
discrete differential operator that calculates the 
approximate value the image brightness gradient. The 
result of applying Sobel operator at each point in the 
image is a vector of brightness gradient at this point, or 
a norm. [ 1 ] 
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Easier to understand, the operator computes the 
gradient of the image brightness at each point. So is the 
direction of maximum increase brightness and the 
change value in this direction. The result shows how 
"sharply" or "smoothly" the image changes brightness 
at each point, and thus the finding a point probability 
on the face and the boundary orientation. In practice, 
the change amount calculation the of of brightness (its 
probability of belonging to a face) reliably and simpler 
to interpret than the direction the calculation. [2] 

Mathematically, the a function gradient of two 
variables for each point in the image (which is a 
brightness function) - a two-dimensional vector whose 
components derivatives are the brightness of the image 
horizontally and vertically. In each point of the image 
gradient vector is oriented in the maximum increase 
brightness direction, and its length corresponds to the 
change in brightness. This means that the Sobel 
operator result at the constant brightness point is the 
zero vector, and at a point on the border different 
brightness levels areas - vector crossing the border in 
the increasing brightness direction. [3] 

Segmentation process using the Sobel operator is 
based on a simple filter mask moving from point to 
point in the image, at each point (x, y) filter response is 
calculated using the originally-formed bonds. In case 
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the linear spatial filter response given by the sum of 
multiplying filter coefficients the corresponding pixel 
value in the region covered by the mask filter. For 3 X 3 
mask element illustrated in Fig. 1 , the result 

R = -I) f(x -\,y-\) + w(-\,0)f(x - 1, y) + ... 

+ w(0,0)f(x,y) + ... + w(1,0)/(jc + \,y) + 

+ w{\,\)f{x + \,y + \) 

(1.1) 

that, apparently, is the sum of the multiplying 
coefficient mask and pixel values directly under the 
mask. In particular, note that the coefficient w (0, 0) is 
at a value of / (x, y), thereby indicating that the mask is 
centered at the point (x, y). Upon detection of 
brightness differences can be discrete analogues of first 
and second order. 

The first derivative dimensional function / (x) is 
defined as the difference between neighboring 
elements values: 

df 



x^The origin of coordinates the image 



dx 



= /0 + i)-/0) 



(1.2) 



Here we used the recording as a partial derivative in 
order to keep the same notation as in the two variables 
fx,y) case, where it is necessary to deal with the partial 
derivatives of the two spatial axes. 

Similarly, the second derivative is defined as the 
difference between adjacent values of the first 
derivative: 



dx 7 



= /(* + !) + /(*-!) -2/0) (1.3) 



The the first derivative calculation of the digital 
image is based on different discrete approximations of 
the two-dimensional gradient. By definition, the image 
gradient / (x, y) at point (x, y) - a vector: 

Gx 



dx 
By 



(1.4) 



As is known from mathematical analysis, the 
gradient vector coincides direction with the maximum 
rate direction of change function f at the point (x, y). 

An important role is has detection the contours this 
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Fig. 1 The spatial filtering scheme, (response) R linear 
filtering point (x, y) images is: 



vector module, denoted V/ and is equal 
V/ = |V/| = yjG 2 x + G 2 Y 



(1.5) 



This value is equal to the change function / 
maximum rate at the point (x, y). 

The gradient vector direction is also an important 
characteristic. Denote a(x, y) the angle between the 
vector Vf at the point (x, y) and the axis x. As is known 
from mathematical analysis 



a(x, y) = arctg 



Gy 
Gx 



(1.6) 



It is easy to find the contour direction in the point (x, 
y), which is perpendicular to the gradient vector at this 
point. And you can calculate the gradient of the image, 
calculating the value of the partial derivatives df /dx 
and df Idy for each point. 

G.v and G Y - these are two of the matrix, where each 
point contains the approximate derivatives in x and y. 
They are calculated by multiplying the matrices G x and 
G y , and the summation of both matrices, as a result of 
the result is written to the current x and y coordinates in 
a new image: 

g = Jg; 

Matrix Gy and Gx: 



f\ + Gy 



(1.7) 
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and 



1 



G, = 



1 



0 0 0 

+1 +2 +1 

-l o +r 

-2 0 +2 
-1 0 +1 



xA 



(1.8) 



xA 



(1.9) 



3. Results and Discussion 

The result of applying Sobel operator is a 
two-dimensional the gradient map for each point. 
This can be process and display a picture, which 
sections to a gradient large value (mainly faces) will 
appear as white lines. [4, 5] The following images 
illustrate this with an segmentation example of medical 
images: 
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Fig. 3. Image using Sobel operator. 
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4. Conclusions 

Thereby demonstrated the practical application of an 
experimental Sobel algorithm to analyze medical 
images. Therefore using Sobel operator in computed 
automated system allows greatly increase diagnostic 
possibilities of medical images segmentation for 
osteoscintigraphy analysis. 

Acknowledgment 

Supported by grant of RFBR N° 13-07-00667 
«Computer-aided analysis of medical images of the 
combined SPECT and X-ray CT» 



References 

[1] K. Engel, M. Hadwiger, J. Kniss, C. Rezk-Salama, 
Real-Time Volume Graphics, A K Peters, Wellesley, MA, 
2006. Pp. 112-114. 

[2] Vatutin E.I., Miroshnichenko S.Yu., Titov V.S. 
Programmnaya optimizatsiya operetora Sobela s 
ispolzovaniem SIMD-rasshireniy processorov semeystva 
x86. In Handbook of Telecommunikatsii, 2006. Ns 6. 
P.12-16. 

[3] B. Jahne, H. Scharr, and S. Korkel. Principles of filter 

design. In Handbook of Computer Vision and 

Applications. Academic Press, 1999. 
[4] Anisimov BV. Raspoznavanie I tsifrovaya obrabotka 

izobrajeniy. - M.: Vyshaya shkola, 1983. 295 p. 
[5] Kudryavtsev L.V. Kratky kurs matematicheskogo analiza. 

-M.:Nauka, 1989. P.736. 



Journal of Mathematics and System Science 4 (2014) 261-268 




PUBLISHING 



On Pairwise Singular Sets and Pairwise Singular Maps 

Rina Pankaj Verma 1 , Mamta Singh 2 

1. School of Studies in Mathematics, Vikram University Ujjain (M.P.), India 

2. Department of Appl. Sci. & Comp. Tech. Bundelkhand University, Jhansi, India 

Received: December 10, 2013 / Accepted: January 16, 2013 / Published: April 25, 2014. 

Abstract: In this paper we introduce the concept of pairwise singular sets and pairwise singular maps between pairwise locally 
compact and pairwise hausdorff spaces and study the properties of pairwise singular maps. 

Keywords: Pairwise singular sets, pairwise singular maps 



1. Introduction 

A Bitopological space is a triple (X, 3i, 32), where 
3i and 32 are topologies on a set X. J. C. Kelly [11] 
initiated the systematic study of such spaces and 
several other authors namely Weston [15], lane [6], 
Patty etc. contributed to the development of the theory. 
Kelly introduced pairwise Hausdorff spaces, pairwise 
regular spaces and pairwise normal spaces in the 
theory [11]. 

Recall that a cover U of (X, 3i, 32) is called 
pairwise open if U c 3i U 32 . and U contains at least 
one non-empty member of 3i and one non empty 
member of 32- If every pairwise open cover of (X, 3i, 
32) has a finite sub cover then the space is called 
pairwise compact [5]. 

According to 1. 1. Reilly [ 9 ] a bitopogical space (X, 
3i, 32) is called pairwise locally compact space if 3i 
is locally compact with respect to 3i and 32 is locally 
compact with respect to 3i. 3i is locally compact with 
respect to 32. if each point of X has a 3i open 
neighbourhood whose 32 . closure is pairwise 
compact. 

A space (X, 3i, 32) is called pairwise Hausdorff if 
for two distinct points x and v there is a 3i 
neighbourhood U of x and a 32 neighbourhood V of y 
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such that U 0 V = c() [11] . 

A function f. (X, 3i, 3 2 ) — > (Y, Xi, X2) is called 
pairwise continuous if the induced functions / : (X, 3i) 
— > (Y, X 1) and/ : (X, 32) — > (Y, X2) are continuous. 

The notion of the singular set of a mapping defined 
by Whyburn and Cain [2, 14] was further investigated 
by various workers including Cain, Chandler, Tzung, 
Magill, Jr. Faulkner and Duda etc.(1981) [2, 3, 5]. 
Different workers described the concept with their 
own point of view, later this concept led to the notion 
of a singular compactification. Infact the 
compactifications which are known as singular 
compactifications are instances of Whyburn's unified 
space. Various authors have seen methods of 
constructing a compactification of X. Given a map f : 
X — > Y where Y is compact authors who should be 
mentioned are Steiner and Steiner, Magill, Blakley, 
Gerlits and Magill, and Chou. Chandler and Tzung 
defined the remainder induced by/ : X — > Y to be L(/) 
= fl {Cly(X - F) I F e K x } and proved that whenever 
Y is compact there is a compactification aX with aX 
- X homeomorphic to L(/"). Cain defined the singular 
set of a map f : X — > Y in the following manner. 
S(/)={peY/VU£ tjip), 
3FeKx,peFcU and/"'(F) £ K x }. 

Chandler and Tzung defined the remainder induced 
by f to be 

L(0=n{Cl Y /(X-F)/F e K x } = {/?eY/VUe 
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77(p),VF £K x> Un/(X-F)^} 
Other definitions of singular set given by Whyburn, 
Chandler, Faulkner etc. of a map between locally 
compact spaces are respectively given as 

S'(f) ={peY/VUe Tjg(P),f-\\j) * K x } 

and 

S(7) ={pe Y/VUe 77k(P), />)^K x }. 

Further it was observed in a paper Singular sets and 
remainders [2] for locally compact and Hausdorff 
spaces that all the definitions taken by various 
workers are equivalent and each of the sets described 
above are actually the same. 

A G-space is a Hausdorff space X on which a 
topological group G acts continuously. The orbit 
space of a G-space X is denoted by X/G and the maps 
which assigns to x in X its orbit is the orbit map on X. 
An equivariant map f from a G-space X to a G-space 

Y induces a continuous map f G on orbit spaces which 
sends the orbit of x in X to the orbit of f(x). Let H be a 
compact subgroup of G and let X be an H-space then 
for h e H and (g, x) e GxX, h{gx) = (gh'\ hx) defines 
an action of H on G X X. The orbit space G x h X of the 
H-space GxX is called the twisted product of G and X 
[13]. 

For H-spaces X and Y, and equivariant map /: X — > 

Y determines a continuous map/^: G x H X — > G ><h Y 
which sends [g, X] in G * H X to [g,fix)] in G x H Y. 

We begin with the study of pairwise singular maps 
and pairwise singular sets and see different definitions 
of pairwise singular set which are equivalent in 
section 2. A brief account of properties of pairwise 
singular maps in relation to product, composition, 
restriction of continuous maps etc. is given in section 
3. The last section 4 is devoted to see the behaviour of 
pairwise singular maps in relation to various structures 
like G-spaces etc. 

For terms and notations not explained here we refer 
the readers to [1, 7, 10, 13]. 

2. Pairwise Singular Sets 

In this section the notion of the singular set of a 



map / between pairwise locally compact, pairwise 
Hausdorff bitopological spaces is introduced. It is 
proved that like topological spaces the four sets 

Sb(0,S b '(0, W) and Z B (f) are the same. 

2. 1 Definition 

Let f. (X, 3i, 3 2 ) -> (Y, Xi, X 2 ) be a pairwise 
continuous map where X, Y are pairwise locally 
compact spaces. Then a point y e Y is called Xi 
singular point with respect to 3 2 if for each open set U 
e Xi of Y contiaining y, 3 2 C1/" (U) is not contained in 
a compact set of X. 

A point p e Y is called a pairwise singular point if 
it is Xi singular point with respect to 3 2 and X 2 
singular point with respect to 3,. 

The set of all pairwise singular points of f is called 
the pairwise singular set of / and it is denoted by 
Sb(/). 

S(f, Xi, 3 2 ) = {p e Y / V U e Xi, 3 2 C1/" 1 (U) is not 
compact} 

S(/,X 2 , 30= (peY/VUeli, S.Cl/^O-O is not 
compact} . 

SbW"= S(f, Xi, 37) n S(/-,X 2 ,30 
Other definitions of singular sets as in case of 
bitopological spaces are the following: 

2.2 Definition 

S'(f, X i, 3 2 ) = {/?eY/VUeX i,/"'(3 2 Cl U) is 
not compact} . 

S'ft Z* 3.) = {p e Y I V U 6 Zi,f U) is 
not compact} . 

s B '(/) = s'(/;Xi,3 2 )n sv,z 2 , 3,) 

2.3 Definition 

S(f, Xi, 3 2 ) = {/?eY/VUe Xi, 3 X 1 compact set 
F with /jgFcU and / _1 (F) is not 3 2 compact}. 

S(/i X 2 , 3,)={£>eY/VUe X 2 , 3 X 2 compact set 
F with p e F cz U and / _1 (F) is not 3i compact}. 

s B (/) = s(/;Xi,3 2 )n s(/;x 2 ,3,). 

and 
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2.4 Definition 

Uf, Zi, 3 2 ) = n { XiCl/(X - F) / F is 3 2 compact}. 
L(f, X 2 , 3.) = (1 { X 2 C1/(X - F) / F is 3i compact}. 

L B (/") = L(/;Xi,3 2 )nx(/;x 2 ,3 1 ). 

First we show that all the sets define above are 
equivalent. 

2.5 Theorem 

Sb(A) = XbOO = S b '(/)= W) 

Proof : If p £ L B (f ), then p * L (f, Z i, 3 2 ) D L (f, 
X 2 , 3,). 

Assume that p g L (f, X i, 3 2 ). 

Then p e n {Xi C1/(X - F) / F is 3 2 compact} 

Then 3 U e Xi such that p e U and U (1 /(X - F) = 
(() for some 3 2 - compact set F. 

Or /"'(U H /(X - F)) = (|) for some 3 2 -compact set 

F; 

Or/- 1 (U)n(X-F) = ((.; 

For any X i compact set K with p e K c U,/"'(K) 
c F it follows that it is compact. 

This implies that p £ S B (f). Thus S B (f) c L b (f) 

Now we show that L B (f ) c S B ' (/")• 

Let p e S' B (A). Then p e S'tf, Zi, 3 2 ) D S'tf Z 2 , 
3,)- 

Assume that /? g S'(/, Xi, 3 2 ), then there exists U e 
Xi with p e U and /"'(3 2 C1 U) is 3 2 compact. 

/(x-/- 1 (u)) = f(x)n(Y- fj). 

peUandUn [f(X) n (Y-U)] =^p ? L B (f, Z u 
3 2 ). 

Next, take p e S B ' CO- Then P e S V, Zu 3 2 ) (1 S'(/", 
X 2 , 3,). 

To show that p e S B (/) , we show that e 



S(^ Z b 3 2 )o 

let p g U, with UeI, choose V with X 2 C1 V c U 
3 2 C1/" 1 (U) (X 2 C1 V). 

If 3 2 C1/" 1 (U) is compact then/" 1 (X 2 Cl V) is also 
compact. 

/? e S(/", Xi, 3 2 ) Similarly it is shown that p 



Hence p e S B (f) implies S B ' (/") c S B (/). 
Finally we show that,/? e S B (f). 

Let /7 e S B (/), take U 6 Xi with p e U. Then 

3 2 C1/ _1 (U) is not compact. 

Choose VcU with p e V 6 X 2 with X 2 CI V c U 

Then 3 2 C1/"' (V) c (X 2 C1 V) 

Note that 3 2 C1 f' l (V) is not compact. This implies 
that/"'(3 2 Cl V) cannot be compact. This implies that 
p e V c I 2 CI V c U satisfying that X 2 CI V is 
compact but f' x (X 2 CI V) is not compact. 

Sop e S B (0- 

This implies that S B (/) a S B (f). 

From above we conclude that 

Sb(A) = Lb(A) = S b ' Cf) = s B (f). 

Throughout this paper, we use the following 
definition of a pairwise singular map. let f: (X, 3i, 3 2 ) 

(Y, Xi, X 2 ) be a pairwise continuous map where X, 
Y are pairwise locally compact, pairwise Haudorff 
spaces. Then a point y e Y is called Xi singular point 
with respect to 3 2 if for each open set U e Li of Y 
contiaining y, 3 2 C1 / _1 (U) is not in a compact set of X. 

A point is pairwise singular point if it is Lingular 
point with respect to 3 2 and X 2 singular point with 
respect to 3,. 

i.e. S B (0 = Stf, Zu 3 2 ) fl S(f, X 2 , 3.) 
where, S(f, Zi, 3 2 ) = {peY/VUe Xi, 3 2 C1/-'(U) 
is not compact} and S(f, X 2 , 3,) = {p e Y / VU e X 2 , 
3,C1 /"'(U) is not compact}. 

2. 6 Example 

Consider the sine map 

sin : (R, u, 3 d ) — > (R, u, 3 d ) 
where w denotes the usual topology on R and 3 d 
denotes the discrete topology on R. The sine map is 
pairwise continuous and since it has non-compact 
fibres, it is a pairwise singular map. 

3. Properties of Singular Maps 

In the present section we discuss a few properties of 
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pairwise singular maps. The behavior of pairwise 
singular maps is seen in relation to the product , 
composition and restriction etc. an analogue of pasting 
lemma is also obtained for pairwise singular maps. 

3.1 Theorem 

Let/: (X, 3i, 3 2 ) -» (Y, Xi, X 2 ) and g: (Y, Xi, X 2 ) 
— > (Z, Xi', X2') be pairwise singular maps. Then the 
composition gof : (X, 3i, 3 2 ) — * (Z, Xi', X2') is a 
pairwise singular map. 

Proof: To show that gof is pairwise singular map, 
we prove the following cases : 

Case I: y e Z is a (Xi 1 , 3 2 ) singular point of gof. 

Case II: y e Z is a (X2', 3,) singular point of gof. 

Case I: First we show that gof is (Xi\ 32) singular. 
Take p e Z and U' e Xi' with p e U'. 

Consider, 3 2 C1 (g o f) A {U) = 3 2 C1 (T V(U'))- 

Since g : Y — > Z is pairwise continuous, g _1 (U') e 
Xi. Since /is pairwise singular, by definition of (Xi', 
32) singular map 32/"'(g" 1 (U')) is not compact in X. 
This implies that gof:X— >Zisa (Xi', 3 2 ) singular 
map. Similarly, it is proved that g o f : X — > Z is (X2', 
3,) singular. 

3.2 Theorm 

Let/: (X, 3 h 3 2 ) (Y, Z h X 2 ) and h: (X, 3,', 3 2 ') 
— > (Y 1 , Xi', X2') be pairwise continuous maps. Then 
the map/x h : (X x X, 3,x 3,', 3 2 x3 2 ') — ► (Y x Y, 
XixXi', X 2 xX 2 ') defined by (fx h) (x, x') = (/(x), h(x% 
where (x, x') 6 X x X' is pairwise singular if and only 
if either f or h is pairwise singular. 

Proof: let / x h : X x X' — > Y x Y' be pairwise 
singular. If f is not pairwise singular then there is a 
point y e Y which is not a pairwise singular point. 
Then there are three cases. 

Case I: y e Y is not a (Xi, 32) singular point of f. 

Case II: v e Y is not a (X2, 30 singular point of f. 

Case III: y e Y is neither (Xi, 3 2 ) nor a (X2, 30 
singular point of f. 

Case I: Since y e Y is not a (Xi, 32) singular point 
of / there exists U e Xi with y e U satisfying that 



SzCl/" 1 (U) is compact. Lety e Y' and U' e Li' with 
y e U'. Since (y, y) e Y x Y is a (X^Xi', 3 2 x3 2 ') 
singular point offxh, we conclude that (32 x 32') CI 
(fx h)~ l (U x U') is not compact. Since (3 2 x 3 2 ') CI (f 
x hf (U x Tj') = 3 2 C1/ _1 (U) x 3 2 ' CI h A (U). 

It follows that 32 1 CI K (U 1 ) is not compact. Thus h 
is a (Xi, 32) singular map. 

The proof of other cases follows similarly. This 
implies that h is a pairwise singular map. The proof of 
the converse is left to the reader. 

3.3 Definition 

Let (X, 3i, 32) be a bitopological space and AcX 
then A is called a pairwise L-Subspace of X if the 
restriction of each pairwise singular map / : (X, 3i, 32) 
— > (Y, Xi, X2) on (A, 3i, 32) is pairwise singular 
where (Y, Xi, X2) is any bitopological space. 

3.4 Theorem 

Let (X, 3i, 32) be any bitopological space, let K be 
3i, 32 compact set of X. Then X - K is a pairwise 
L-subspace of X. 

Proof: let f : X —> Y be a pairwise singular map, 
where Y is pairwise locally compact space. 

The restriction of / : (X, 3i, 3 2 ) — > (Y, Xi, X2) is 
denoted by g : (X - K, 3 b 3 2 ) (Y, Xi, X 2 ). To 
prove that g is a pairwise singular map we show that 

Case I: g is a (Xi, 32) singular map. 

Case II: g is a (X2, 3i) singular map. 

To show that g is a (Xi, 32) singular map, take a 
non empty open set U e Xi 

Case I: If g "'(U) =/ _1 (U), then since /is pairwise 
singular map, g"'(U) is not contained in any 32 
compact set of X - K. 

Case II : If g A (U) ^f'\U) then/ _1 (U) c g\\J) u 

K. 

If g~ (U) is not a (Xi, 3^) singular map then g~ (U) 
is contained in a 32 compact set F of X - K. 

Then/ _1 (U) c F u K which is a union of two 32 
compact sets of X. This implies that / "'(U) is 
contained in a 32 compact set of X which is a 
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contradiction. Hence g is a (Zi, 3i) singular map. The 
proof of case II is similar. 

Since the pairwise singular map and the 
bitopological space Y are arbitrarry, X - K is a 
pairwise L-subspace of X. 

3.5 Theorm 

Let X be a bilopological space and A, B be two 
pairwise L-subspace of X with either A or B 3 1, 32 
closed in X. Then A u B is a pairwise L-subspace of 
X. 

Proof: let f : (X, 3 b 3 2 ) -» (Y, Zi, Z 2 ) be a 
pairwise singular map, where Y is a pairwise locally 
compact bitopological space. Assume that A is a 3 1, 
3 2 closed L-subspace of X. f A , fa and f AuB denote 
respectively the restrictions of f on A, B and A uB. 
Since A and B are pairwise L - subspaces, f A and f*B 
are pairwise singular maps. K be a 3 1, 32 compact set 
of A uB. To prove that f AU B is a pairwise singular 
map i.e. f AuB : (AuB, 3 b 3 2 ) -> (Y, Zi, Z 2 ) is a (Xi, 
32) and (X 2 , 3;) singular map, we proceed as below. 

First we show that f AwB is (Xi, 32) singular map. 
Take U e Zi. Then /" Au b(U) = /"a(U) u /"' b (U) 

If /^aubCU) is contained in a compact set F of 32 in 
AuB then / _1 A (U) is contained in F fl A. Since F fl 
A is a closed subspace of F, it is compact. Since / 
_1 A (U) c F fl A, we have a contradiction to the 
assumption that f A is a pairwise singular map. 
Similarly it is proved that/vuB : (A u B, 3i*, 3 2 *) 
— > (Y, Xi, X2) is a (X2, 3i) singular map. Since the 
map / and the bitopological space Y are arbitrary, A u 
B is a pairwise L-subspace of X. 

3. 6 Theorm 

Let X be a bitopological space and let A be a 3 1, 32 
closed L-subspace of X. If B is a subset of X such that 
AcB, then B is also an pairwise L-subspace of X. 

Proof: let f : X — > Y be a pairwise singular map 
and f A and f B be restrictions of / on A and B 
respectively. Since A is pairwise L-subspace of X,f A 
is a pairwise singular map. 



To prove that / B is a pairwise singular map, we 
show that f B : (B, 3j, 32) — > (Y, Xi, X2) is pairwise 
singular map i.e. it is a (Zi, 32) and (X 2 , 3i) singular 
map. First we prove that f B is (Xi, 32) singular map. 
Take an open set U e Xi in Y. Since A cB, 
f-\(\J) c/-' B (U). 

If J~ 1 b(U) is contained in a compact set F of B then / 
_1 A (U) is contained in F fl A of 32 which is a closed 
set of a compact space and hence compact. This 
contradicts the hypothesis that/X is a pairwise singular 
map. Hence f B is a (Xi, 3 2 ) singular map. Similarly it 
is a (X2, 3i) singular map. 

Since the pairwise singular map / and the space Y 
are arbitrary, B is a pairwise L-subspace of X. 

The following proposition gives an analogue of 
Pasting lemma for pairwise singular maps. 

3. 7 Theorm 

Let X, Y be pairwise locally compact spaces and A, 
B be 3i, 32 closed sets of X such that AuB = X. let f : 
(A, 3!', 3 2 ') - (Y, Xi, X 2 ) and g : (B, 3r, 3 2 ") - (Y, 
Zi, Z 2 ) be two pairwise singular maps such that/[x) = 
g(x) forx 6 AflB. Then the map h : (X, 3i, 32) — » (Y, 
Zi, Z2) defined by 

w Xg(x):ifxeB 
is pairwise singular. 

Proof : To prove that h is pairwise singular, we 
show that h is (Zi, 32) singular and (Z2, 3i) singular. 
First we show that h is (Zi, 32) singular, let U be a Zi 
open set of Y. 

Then h-\\J)=f- l (U)ug 1 (y) 

or 

3 2 C1 h'\\J) = 3 2 'Clf l (U) u 3 2 "C1 g 4 (U) 
Since f and g are pairwise singular maps, and 32'Cl 
/"'(U) and 32"Clg _1 (U) both are non compact. Hence 
32Cl/z" 1 (U) is also non compact. Since the open set U 
e Zi in Y is arbitrary, h is a (Zi, 32) singular map. 
The proof of case (ii) i.e. h is (Z2, 3i) singular map 
follows similarly. This implies that h is a pairwise 
singular map. 
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4. Pairwise Singular Maps And Various 
Structures 

In this section the behaviour of pairwise singular 
maps is seen in relation to the G-spaces, Cone and 
twisted product etc. 

4. 1 Definition 

A pairwise G-space consists of the following : 

A Bitopological space (X, 3i, 32), a bitopological 

group (G, Zi, Z 2 ) and a pairwise continuous map 9: G 

xX->X satisfying 
9 (e, x)=x 

9 (gu 9fe, x)) = Q(gig 2 , x) V x e X, g u g 2 e G. 
As in the case of a single topological space, the 
point 9(g, x) is denoted by gx. The pairwise orbit of a 
point x 6 X is denoted by Gb(x) = {gx I g e G} the 
collection X/G of pairwise orbits is given two 
topologies and thus becomes a bitopological space. It 
is called the Bi-orbit space of x e X. The map q : X 
— » X/G taking x to its orbit is called the pairwise orbit 
map. 

4.2 Theorm 

Let X, Y be pairwise G-spaces with G-3i, 32 
compact bitopological group and / : X -> Y be a 
pairwise singular equivariant map. Then the induced 
map fa : X/G — » Y/G is pairwise singular. 

Proof: Consider the following commutative 
diagram : 

/ 

(X, 3„ 3 2 ) ►(Y, Zi, X 2 ) 

Ixi i q y 

(X/G, 3/, 3 2 ') ►(Y/G, Zi', Z 2 ') 

/g 

where <7. v and q y are pairwise orbit maps. 

To show that /q is a pairwise singular, we prove the 
following 

Case I : G(y) e Y/G is a (Zi', 32*) singular point of 

fa. 

Case II : G(y) e Y/G is a (Z2', 3i') is a singular 
point of fa. 



Case I : Suppose G(y) e Y/G is not a (Zi', 3 2 ') 
singular point of fa- Then there exists an open set U of 
Zi' containing G(y) such that 32'Cl fa~ X (\l) is compact. 

Since q x is a 32'-compact map, ^ a " 1 (3 2 'C1(/g" 1 (U))) 
is 32-compact. From continuity of q x and 
commutativity of the above diagram it follows that 

q ;\^c\(f G \\j))) => 3 2 ci 9 ; 1 (/b" 1 (u)) = SaCic/" 1 ^ 

"'(U))) 

Thus 32Cl(/" 1 (g' r " 1 (U))) is 32 compact which is a 
contradiction to the hypothesis that f is (Zi, 32) 
singular. 

The proof of case II follows similarly. This implies 
that fo is a pairwise singular map. 

let (G, 3i, 3 2 ) be a bitopological group. Assume 
that G is 3 t , 3 2 Hausdorff. let H be a 3 19 3 2 compact 
subgroup of G and (X, Zi, Z2) be a pairwise H-space. 
Then for h e H and (g, x) e G x X, h(g, x) = (gh~ l , hx) 
defines an action of H on G x X which is pairwise 
continuous. The pairwise orbit space (G x H) / H is 
called the pairwise twisted product of G and X. For 
pairwise H-spaces X and Y, a pairwise equivariant 
map / : (X, Zi, Z2) -» (Y, Zi', Z2') determines a 
pairwise continuous map f ti : Gx H X -> Gx H Y which 
sends [g, x] e Gx H X to \g,fix)\ in G X H Y. 

4.3 Corollary 

Let f : X — > Y be a pairwise singular map then^H : 
G X H X — > G X H Y is also a pairwise singular map. 

Proof: Consider the following commutative 
diagram 

I G x/ 

(G x X, 3 l5 3 2 ) ►(& x Y, Zi, Z 2 ) 

q x 4 I q y 

(G x HX, 3/, 3 2 ') ►( G x HY, Zi', Z 2 ') 

h 

where q x and ^ ; are pairwise orbit maps. 

We have to show that f H is a pairwise singular map. 
i.e. 

Case I :f H is (Zi', 3 2 ') singular. 

Case II :f H is (Z2', 3i') singular. 

Case I : Suppose f H is not (Zi', 32 1 ) singular. Then 
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there is an open set U e Xi' in G x H Y such that 3^0 
/^"'(U) is compact. Since q x is a 3 2 ' compact map, 
q x A (^2ClfH '(U)) is also compact. From continuity of 
q x and commutativity of the above diagram, it follows 
that 

q;\3 2 'C\(f H -\U))) 3 3 2 Cl 9 ; 1 (fe" 1 (U)) = W(qy 
0(IG xf)f(U). 

Since 32C1(<7,,. o (Ig x f)) " (U) is 32 compact, a 
contradition to the hypothesis that (Iq x f) is pairwise 
singular is obtained. This implies that f H is (Xi', 32 1 ) 
singular. The proof of case (II) follows similarly. Thus 
fn is pairwise singular map. 

4.4 Definition 

Let (X, 3i, 3 2 ) be a bitopological space. The 
pariwise cone (TX, Xi, X2) is the quotient 
bitopological space X x I/A, where A = X x { 1 } . 

4.5 Theorem 

Let f : X — > Y be a pairwise singular map where X 
and Y are pairwise locally compact spaces. If the 
quotient map p : X x I — > TX is a 3i, 32 compact map, 
then the induced map Tf : T x -> T Y is a pairwise 
singular map. 

Proof: Consider the following commutative 
diagram 

(X x I, 31, 3 2 ) ► (Y x I, Xi, X 2 ) 

pi iq 

(TX, 3i\ 3 2 ') ► (TY, Xi', X 2 ') 

T/ 

Now, we have to show that Tf is pairwise singular, 
i.e. 

Case I : T f is (Xi 1 , 32 1 ) singular. 

Case II : T f is (X2', 3i') singular. 

Case I : Suppose that Ty is not (Xi', 32*) singular. 
Then there exists an open set U e Xi' of Ty satisfying 
that 32 , C1T/" 1 (U) is compact. Since p is a 32' compact 
mapping jp" 1 (32 i C1T/" 1 (U)) is 32-compact. From 
continuity of p and commutativity of the above 
diagram, it follows that p l (3 2 'ClTf\\J)) 3 3 2 C1 



(p-\T f - l (U))) = 3 2 C1 [{fx l)-\q- l (V))] Thus 3 2 C1 (fx 
I)" q" (U) is compact, a contradiction to the hypothesis 
that/is a pairwise singular map. This implies that T/is 
a (Xi', 32') singular. The proof of case (II) is similar. 
So Tf is a pairwise singular map. 

5. Discussion and Conclusion 

We discussed properties of pairwise singular maps. 
The behavior of pairwise singular maps in relation to 
the product , composition and restriction etc. The set 
of all pairwise singular maps from R to R which are 
periodic in addition is closed under the 
addition,product ,division and composition of maps. 

References 

[I] G. E. Bredonjntroduction to compact transformation 
groups, academic press, Newyork (1972) 

[2] George L. Cain, Richard E. Chandler And Gary D. 
Faulkner Singular Sets and Remainders, Trans. Amer. 
Math. Soc. 268, 161-171. 

[3] Richard E. Chandler, Gary D. Faulkner, Joshephine P., 
Guglielmi and Margaret, Memory Generalizing the 
Alexandroff Uryshon Double Circumference 
Construction, Proc. Amer. Math. Soc. 83 (1981), 
606-608. 

[4] Richard E. Chandler and Gary D. Faulkner, Singular 

Compactification: The Order Structure, Amec. Math. Soc. 

100(1987), 377-381. 
[5] R.E.Chandler, Hausdroffcompctification, Marcel Dekker, 

New York (1979). 
[6] J. Dugundji, Topology Allyn and Bacon, Boston (1966). 
[7] Gary D. Faulkner, Compactifications whose remainders 

are retracts, Proc. Amer. Math. Soc. 103 (1988), 

984-989. 

[8] Josephine P. Gaglielmi, Compactifications with Singular 
Remainders, Ph.D. Thisis North Carolina State 
University. 

[9] J.C. Kelly, Bitopological spacecs, Proc London Math. 

Soc 13(1963), 71-89. 
[10] Y.W. Kim, Pairwise compactness, Publ Math, Debrecen 

15(1968), ,87-90. 

[II] E.P.Lane, Bitopological spaces and Quasi-uniform spaces, 
Proc.London Math.Soc Vol.17 (1967), 241-256. 

[12] K.D. Magill, Jr., A note on compactifications, Math. Z. 

9(1966) 322-325. 
[13] James R. Munkres, Topology: A First Course, Prentice 

Hall of India Pvt. Ltd. (1983), New Delhi. 
[14] C.W. Patty, Bitopological spaces, Duke Math, J. Vol. 



268 



On Pairwise Singular Sets and Pairwise Singular Maps 



34(1967), 387-392. 
[15] I.L. Reilly, Bitopological Local Compactness, Indag. 

Math 34 No 5(1972), 407-410. 
[16] I.L. Reilly On Pairwise connected Bitopological spaces, 

Kyungpook Math J. 1 1(1971), 22-25. 
[17] Anjali Shrivastava and Abha Khadke, On Singular and 

H-Singular Maps, Bull. Cal. Math. Soc. 90(1988), 

331-336. 

[18] P. Srivastava and K.K. Azad Topology Vol I, Shivendra 



Prakashan, Allahabad (1985). 
[19] J.D. Weston, On the comparison of Topologies, J. 

London Math. Soc. 32(1957), 342-354. 
[20] G.T. Whyburn, A Unified Space for Mappings, Trans. 

Amer. Math Soc. 74,844-850. 
[21] G.T. Whyburn, Compactification of mappings, Math. 

Ann. 166(1966), 168-174. 
[22] Eric K. van Douwen, Retraction from fSX onto |3X\X, 

General Topology Appl. 9(1978), 169-173 



Journal of Mathematics and System Science 4 (2014) 269-274 



The Convergence of 1 -Periodic Branched Continued 
Fraction of the Special Form in Parabolic Regions 
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Abstract: Branched continued fractions are one of the multidimensional generalization of the continued fractions. Branched continued 
fractions with not equivalent variables are an analog of the regular C-fractions for multiple power series. We consider 1 -periodic 
branched continued fraction of the special form which is an analog fraction with not equivalent variables if the values of that variables 
are fixed. We establish an analog of the parabola theorem for that fraction and estimate truncation error bounds for that fractions at 
some restrictions. We also propose to use weight coefficients for obtaining different parabolic regions for the same fraction without any 
additional restriction for first element. 

Key words: Continued fractions, 1 -periodic branched continued fraction of special form, convergence, uniform convergence, 
truncation error bounds. 



1. Introduction 

Establishing the convergence regions is one of the 
investigation methods of convergence for continued 
fraction 

a„ . a, 



l + D- 

n=l 1 



1 + - 



i a 2 

1 + — 

1 + . 



(1) 



where a n e C . The sequence of nonempty sets 
{Q n }^ =] , Q n c C are called the sequence of 
convergence regions for fraction (1), if conditions: 
a n 6 Q. n for n > 1 guarantee the convergence of that 
fraction. If Q n = Q ( n > 1 ), then the region Q 
called the simple convergence region of continued 
fraction (1). It is often considered the circular, angular, 
parabolic regions in the analytic theory of continued 
fraction. 

In particular, the research review concerning the 
parabolic regions is given in the next monographs [10, 
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13, 14, 15]. 

We consider multidimensional generalization of 
parabolic theorem for 1 -periodic branched continued 
fraction (BCF) of the special form [10, p. 99]. 

Analogs of parabolic theorems are established for 
different forms of branched continued fraction 
(multidimensional generalization of continued 
fraction). D.I. Bodnar [4] proved multidimensional 
analogs of parabolic theorems for BCF of general 
form with iV branches 



00 N a : 



1+Z> E 



m 



k=\ 



h =i 



N 

i 1= l 1+ £ a i{2) 



(2) 



where a i(k) e C , i(k) = i l i 2 ...i k - multi index 
(l<i k <N, k>l)- T- M - Antonova [1] established 
multidimensional generalization of parabolic theorem 
with the following condition: if elements from 
different "floors" belong to different parabolic regions. 
R.I. Dmytryshyn [9] proved the convergence of 
multidimensional generalization g-fractions in the 
union of parabolic regions and established the 
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truncation error bounds for that fractions in some 
bounded part in the parabolic region. O.Ye. Baran [3] 
investigated the convergence of branched continued 
fraction of the special form 



* =1 i t =l»i(*) «i=l 



-, (3) 



'^(D+Z 

,2=l6 «2)+.. - 

where a i(k), fy(k) - complex numbers, i(k) - multi 
index, \<i k <i k _ { , k > 1 , i 0 = N - integer, if 
partial numerators are in corresponding parabolic 
regions and partial denominators - the half-plane. 
Kh.Yo. Kuchmins'ka [11] proved analogs of parabolic 
regions for two-dimensional continued fractions with 
complex elements 

00 a i i °° a i+i i 00 a i i+ i 

k=i®, j=i 1 j=\ 1 

2. Main Definitions and Methods 

The object of our researching is 1 -periodic branched 
continued fraction of special form, that we receive 
from (3), if we put: a i(k) =c (j) b i(k) =1 (l<i k <i k _ x , 
k > 1 ), that is BCF next form 

f V 1 



CO h-1 C: ^ 

1+Dlf 

k=M.=\ 1 



i+Z- 



'l+z 



(4) 



•J 



; 2 =il + 

where Cj - complex numbers for j = l,N , i 0 = N - 
integer. The re-fh approximant of 1 -periodic BCF (4) is 



n <i-l C: ^ 



(n>l; F 0 =l). 



We define 

, . n h-\ C : Jo 

*=Ui=l 1 71=1 



i+ z 

72 = 



1+- fe 1 c >„ 



as n-th tail g-fh order of 1 -periodic BCF (4) (q = l,N ; 
n > 1 ; 7o = ? : fit) 9 ' = 1 ; ^ = 1 )• Obviously, that 



the tails R[ q) ( n >\ \ q = 2,N ) satisfy following inequalities 



recurrence expression 



(?) 



R 



(i) 



(5) 



We use 

• technique of value set for establishing the value 
region of fraction (4) [12, p. 110]; 

• multidimensional analog of Stiltjes-Vitali 
Theorem from [4, theorem 2.17, p. 66] for proving 
uniform convergence of 1 -periodic BCF; 

• method of fundamental inequalities for 
establishing truncation error bounds [15]. 

3. Main Results 

3.1 Useful Lemma 

Let assume that the vector v = (v l ,v 2 ,...,v N ) is 
given, where v ; > 0 - weight coefficients (j = l,N) 

and Z7=i v y < 1 - We denote &j = 1 - TL v s .Then 

following lemma is valid. 

Lemma 1. Let elements Cj ( j = 1, N) of fraction (4) 

satisfy: c, e Pj if) for y e (- f ; f ), where 

P^OO = |z g C : I z I -9*(ze -2!> )< 2p ■ cos 2 (6) 
and p j =v j 3 j . Then tails fi^of 1 -periodic BCF 

O>0,j=l,/V) belong to the half-planes # (.9; J/), 
where 

H(#j;y)= {zeC: ^(ze~ ir )> 9 j cos/}. (7) 

Proof. We consider the half-planes 
\ + H{-v s ;y) = H{& j ;y) for any j = IN . Since 

2j;=i v i ,< 1 next conditions are valid: 0e#(,9,-;r) 
The set Cj/H^y) belongs to the circle 



Q, 



Including 



zeC: 



-'T 



29j cos y 



23 j cos y 



Q. j {y)^H\-v j ;y) i s holding, if 



2i9,'Cos y 



< Vj COS ^ + SR 



fl ( cos y 



are valid. 
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They are equivalent to c ; e Pj (y) ( j = 1, N ). Q.E.D. 
3.2 Main Theorem and Consequence 

We denote G = [z e C : |arg(z + 1]| < 

Theorem. Let elements c ; of fraction (4) belong to 

parabolic regions: P/y) ( j ! = ^>N), which are defined 

by formula (6). Then 

(1) 1 -periodic BCF (4) converges uniformly on any 
compact of the set P, where P = P l (y)x...xP N (y); 

(2) the value region is the circle 



D = \z e C N : |z -| < r] on me set Obviously, that 

D a P . The functional fraction converges on the 
compact D by the multidimensional analog 
Worpitzky's Theorem [2]. It proves the uniform 
convergence of 1 -periodic BCF (4) on any compact of 
the set P. 

We use the following inequality for the truncation 
error bounds of (4) 



\F - F \< 

n+m n \ — 



l 



I 1*1 I 1*2 

c, C, " 



R (N J R (N) \ * 1+ -.+T„=„+i 



n+m n 



k,>o(j=\,N) nn 

7=1 r=l 



R 



U) 



R 



U) 



K N (y) = \ Z £C: 



29 N cosy 



; u v» i 5(9) 'Z"^ 0 

2^ cosy}' where *,-=»-I/=,- + i*/ ; = J J B = _ 1 - 



(3) if \cj\<3j 2 as 2 r for j = 2,N, then the 
truncation error bounds of (4) hold 



\p - F < J 

\ r n+m r n — ^ 



N-l } 
n + N-1 



p n+l (n>0,m>0), 



if 

that leads from difference between F n+m and F n (n>0, 
m>0). It was proved in [7]. 

Since tails is ra-th approximant of 1 -periodic 
continued fraction, c x e P v (y) and P y (y)czG , we 



obtain R ( „ 



for n > 0 



where 



where L = — r — — , p = maxlp ■ ) , p, 

j=l,N " 



ffj 2 cos 2 r 



l-Jl+4c, 



l+./l+4c, 



1c j 1 ,/ 2 1+ V 1+4e i I 

p ; = 2 ' , , Mi =- L - : L , M,=max)l; — ^ — . 

Proof. We use multidimensional analog of 
Stiltjes-Vitali Theorem for proving uniform 
convergence of 1 -periodic BCF. We are going to 
investigate the functional fraction following form 

1+Dlf 

and it's respective the sequence re-th approximants 

{F„(z)}™ =l mdz=(zi,Z2,—, z N )- 

According to the previous lemma it holds 

Ri N) (z)^H(& N ;y). Since F n (z) = (R ( n N \z))~\ then 
F n (z)e K N (y) , that is why the sequence of 
approximants F n (z) is bounded uniformly. We define 

R= rcm{^ ;Vj9j} and the compact D, where 

j=l,N 



x = 1(1 + 71 + 4^") , y = |(l-^+4q) - the 
attracting and the repelling fixed points of the linear 
fraction transformation s{a>) = \ + c l l a> . We use for 
proving the same scheme as in problem 13 [13, p. 49]. 
Otherwise, it holds c l =x- y . Since 

n^K^^f 1 l{\-\ylx\) 

for or q = Sj thus we obtain for 

1 < k x < n + 1 



( 1 1 A 



-1 



n H* (1) 

A X \j Sj +m— r I I ,¥j — r 

where 



<M, 



i2 

vi x l y 



<M, 



zl 

vi x ly 



<M lP ? 



21+Jl+4c, 

M 1= ^ 

1 1-p, 



Pi = 



l-Jl+4c, 



, VT = i. 



l+^/l+4e, 

According to the lemma the next estimations hold 
for modulus of the tails: \Rl n \ > ^{r^ e~ ir )> 3j cosy. 
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Using the inequalities: 
obtain 









-r\\ & lj-r\ 







we 



k+l 

n 



where 



<M jP ) +l {\<k<n + \), 



R 



U) 



R 



U) 



M;=maxjl; — ^ — , p.=— ^ 



cos y 



We denote p = max \p ■ } and obtain the final 
estimations of (4) 



j=l,N 



\ F n+m F "\- \ R (N)\\ R (N)\ 



I Mj/jj ...M N p N N <L 

k i +---+k N =n+\ 



N-l 
n + N-l 



n+1 



where L = U» M Mjfyi™\.^. Q.E.D. 

This theorem generalizes the theorem 2 [6] if we 
define v in such way v = {jfr, jfr)- 

If element q e G and it is fixed this theorem can 
be improved. 

Consequence. Let N>1 and elements Cj of the 
fraction (4) satisfy the conditions: q e G , 

Cj g Pj («! ) , where 



where p. = ^yfe- and the angle ai is 

argq if argq * ;r 



2(AT-1) 

defined as 2a, = 



, £ - some 



0, if axgc l =n 
parameter (0<s<l). Then 

(1) 1 -periodic BCF converges uniformly for c 2 , 
c w if element c } is fixed; 

(2) the value region is the circle 



K = \ Z eC: 



z-- 



£cosa t 
Proof. We take v = (- 1 



£cosa 



i j 



| . Since 



V2 ' 2(N-l) ' ' 2{N-\)I 

qeG and G = (J re (_| ; |)^i<7) for £>7=i/4, then 



3aj : q e ^(a^ . We choose and fix the element cj. 
Since the main theorem with the conditions: 
Cj 6 Pj (a x ) for j = l,N and obtain, that fraction (4) 
converges uniformly on any compact of the set P. 
We define R = max |c and the compact D of the set 

P 

D = {z = ( Zl ,...,z N )eC N : Zj e P.(a t )n \ Zj \< r}\ 

The fraction (4) converges on the compact D. The 
uniform convergence of BCF (5) according c 2 , c N 
leads from the multidimensional generalization of 
Stiltjes-Vitali Theorem. 

Taking into consideration, that 

9 N = 1 - 2^=2 Vj=s/2,we obtain the value region is 
K(s,ai). It is following from the main theorem. Q.E.D. 

It was proved in [5] that F = (fj J =1 Xj ) , where Xj - 
the attracting fixed point of the linear fractional 



transformation of s , = 1 + 



( j = 2,N ) and 



x l = ^ (l + *Jl + 4q ) - the attracting fixed point 
respective of q = 1 + q / co . It is equivalent 
F=l/u N , where - the attracting fixed point of 
t n (co) = u N _ x +c N I co andw;=x/. 



Pj{a x ) = {zeC:\z\ -^(ze^ )< 2 Pj cos 2 «,} 3J fiwonpfe* 



Example 1. Let jV=2, then we consider 1 -periodic 
branched continued fractions of the following form 



1 + 



- + - 



1 + 



1 + ^L 



1 + 1 + - 

1 + — — 1 + — — + - 2 



1 + 



1+ 



1+ 



1 + 



where cj=l+i and c 2 =i/8. The convergence of this 
fraction leads from the main theorem for y=n/4 and 

v = (l/2;l/4) . The calculation is in the next table 
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n 


R d) 


R <2) 




0 


1 


1 


1 




14 


1. 69388+0.41847/ 


1.71306+0.48626; 


0.54022-0. 15335/ 




19 


1.69389+0.41849/ 


1.71303+0.48632; 


0.54021-0.15336;- 



We are going to calculate the attracting fixed points: 
u h u 2 and the value F. It was proved in [5], that 

lim = u x , where u x = ^ (l + + 4q ). Otherwise, 



uj is the solution of the equation u[ 
attractive fixed point of 
transformation r 2 (&>) 



u x +c 2 l Q) 



-u l -c l = 0 . The 
linear-fractional 
is equal 



V"i 



+ 4c 2 ] . That is why the value of 



BCF with 2 branches equals 

F=(u 2 y 



— + — Jl + 4c, +—J— + C, +4c, + — A /l + 4c, 
4 4 V 1 2\2 1 2 2^ 1 



And now we obtain: U!=1.6938+0.4187i for n=14 and 
u 2 =1.71303+0.4863i for w=i9. Towards the end, the 
value of that fraction is equal F=0.5402-0.15336i. 

Example 2. Let a=i/6 and c 2 =10+lli. The 
convergence of that fraction leads from the main 
theorem for y=n/4 and v = (l/2 ;l/4) . 

Analogically as in the previous example we obtain 
Ul =l. 0246+0. 1588i for n=6 and u 2 =4.0728+1.6355i 
for n=58. The value of the fraction equals 
F=0.2114-0.0849i. 

4. Conclusions 

The uniform convergence of the 1 -periodic 
branched continued fraction of the special form is 
proved in the parabolic regions. Truncation error 
bounds are established at some additional restrictions 
of the elements beginning from the second. We 
obtained the value set and the value of the 1 -periodic 
branched continued fraction of the special form. 
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Abstract: In this paper, we consider Variational Iteration Method (VIM) and q-Homotopy Analysis Method (q-HAM) to solve the 
partial differential equation resulted from Fingero Imbibition phenomena in double phase flow through porous media. We further 
compare the results obtained here with the solution obtained in [12] using Adomian Decomposition Method. Numerical results are 
obtained, using Mathematica 9, to show the effectiveness of these methods to our choice of problem especially for suitable values of h 
and n. 
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1. Introduction 

We consider the following nonlinear partial 
differential equation resulting from Fingero-Imbibition 
phenomena in double phase flow through porous media 
[12] 

> 2 Q2« 



as 

dt 



{dx) 



+ S 



d l S 
dx 2 



(1) 



with appropriate initial condition. 

Many years ago, various methods have been put to 
use successfully to obtain analytical solutions of both 
ordinary differential equations and partial differential 
equations such methods as Adomian Decomposition 
Method (ADM) [12], Variational IterationMethod 
(VIM) [1, 8, 20], Homotopy Perturbation Method 
(HPM) [7], and EXP-function Method [19]. One of the 
effective analytical approach to solving nonlinear 
differential equations is Homotopy Analysis Method 
(HAM) [13, 14]. However, the region of h of 
convergence is still questionable as to how large this 
region could be for fast convergence. 

This paper therefore considers a more general form 
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of HAM called q-HAM [6,17] and VIM, for effective 
comparison, to solve equation (1). Recently, q-HAM 
was succefully applied to obtain an aproximate 
solution of Riccati equation in [6] and it was proven 
that the convergence region of series solutions obtained 
by q-HAM is increasing as q is decreasing see also [17, 
18]. Analytical solution of equation 1 is obtained using 
VIM and q-HAM subject to some appropriate initial 
condition. We further compare these results with the 
result obtained using ADM in [12] to affirm the 
reliability of these methods including numerical values. 
We obtain numerical results at the end with different 
values of ft and h ( ft is the new parameter 
embedded in q-HAM) given the effect of both 
parameters on solution of fingero-imbibition 
phenomena in double phase flow through porous 
media. 

2. The Idea of Variational Iteration Method 

We present the basic idea of the Variational Iterative 
Method (VIM) in this section by considering the 
following nonlinear partial differential equation 

d"'u(x,t) 



dt" 



■ + R[u(x,t)] + N[u(x,t)] = f(x,t) (2) 



where R is a linear operator, is a nonlinear 
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operator and f(x, t) is a given continuous function. 

The idea is now to construct the following correction 
function for (2) 



U n+l = U n + 



-f Mt,T) 



d m u(x,t) 
dt m 



+ R[u(x,t)] + N[u(x,t)]- f{x,t) 



(3) 

where A — A(t,z) is called the Lagrange multiplier 
which can be identified optimally by variational 
calculus and u n is the nth term approximate solution. 

The Lagrange multiplier X in (3) of equation (2) ia 
calculated using 

, (-\) m {T-ty- x 

A = - — — — . (4) 

(m-1)! 

See [9, 10] for the interpretation and determination 
of various Lagrange multipliers. 

3. The Idea of q-Homotopy Analysis Method 

Considering the following differential equation of 
the form 

N[u(x,t)]-f(x,t) = 0 (5) 

where iV is a nonlinear operator, (x,t) are 
independent variables, f(x,t) is a known function 
and u(x,t) is an unknown function. 

(l-nq)L(®(x,t;q)-u 0 (x,t)) 

= qhH(x,t)(N[0(x,t;q)]- f(x,t)), (6) 

where n > 1 , Q e [0, 7M denotes the so-called 
embedded parameter, L is an auxiliary linear operator, 
h ^ 0 is an auxiliary parameter, H{x,t) is a non-zero 
auxiliary function. 

Clearly, when q = 0 and q = \ , equation (6) 
becomes 

<S(x,t;0) = u Q (x,t) and ®(x,t;—) = u(x,t)(i) 

n 

respectively. 

So, as <1 increases from 0 to 7, the solution 
0(x,t;q) varies from the initial guess u 0 (x,t) to the 



solution u(x,t). 

If u 0 (x,t) , L , h , H(x,t) are chosen 
appropriately, solution (^{x,t',q) of equation(6) 
exists for q e [0, ^-1 . 

Taylor series expansion of 0(x,t;q) gives 



$>(x, t; q) = u 0 (x, 0 + X M m ^ W 



(8) 



m=l 



where 



u m (x,t) = 



1 d m O(x,t;q) 



ml 8q" 



(9) 



q=0 



We suppose that the auxiliary linear operator L, the 
initial guess u Q , the auxiliary parameter h and 
H(x,t) are properly chosen such that the series 8 
converges at q — \, then we have 



(10) 



u(x, t) = u 0 (x, t) + £ u m (x, t) 

m-1 

Define 

u n ={u Q (x,t),u l (x,t),---,u n (x,t)}. (11) 

Differentiating equation (6) m -times with respect 
to the (embedding) parameter q , then evaluating at 
q = 0 and finally dividing them by m\ , we have the so 
called mth-order deformation equation see Lioa [13, 14, 
16] as 

L \u m (x, t) - X* m u m _ x (x, 0] = hH(x, t)R m (u m -i) • 

(12) 

with initial conditions 

M «(x,0) = 0, k = 0,1,2,..., m-1. (13) 

where 

- , u 1 d m - l (N[®(x,t;q)]-f(x,t)) 
R m {u m -i) ^ m _ 1 | 9=0 

(14) 

and 

. f0 m<l 

Z m =\ fU . (15) 

[n otherwise 

Remark 1. It should be emphasized that U m (x,t) 
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for m > 1 , is governed by the linear operator (12) with 4.2 Application of q-HAM to Fingero-Imbibition 
the linear boundary conditions that come from the phenomena 



original problem. The existence of the factor (^j" 

gives more chances for better convergence, faster that 
the solution obtained by the standard HAM. Ofcourse, 
when n = 1 , we are in the case of the standard HAM. 

4. Applications 

4.1 Application of VIM to Fingero-Imbibition phenomena 

We consider equation (1) with the initial condition 

S(x,0) = \-%e x 

where £=1.11. This problem was solved in [12] 
using ADM. 

In the case of VIM, we use 
S 0 (x,t) = S(x,0) = l-^e~ x as the initial 
approximation and 

m = l, R[S(x,t)] = 0, 



N[S(x,t)] 



ydxj 



d 2 S 

dx 2 



and f(x,t) = 0. 



Hence, using (3), we have the correction function as 



U n+ l =U n +j' Q A(t,T) 



dS(x,t) 
dz 



ydxj 



d 2 S 
dx 



dr 



(16) 



By (4), we obtain A — — 1 . Therefore, we get the 
folio wings from the iteration formular given by (16) 
5, (x, 0 = 1- %e~ x + 2% 2 e- 2x t - %e~ x t (17) 

S 2 (x,t) = 



1 



x.2 



-9£V 3 V -6£V 3 ¥ + ^ 2 e 2x t 3 + W. 



(18) 



Similarly, S m (x, t) for m = 3, 4, 5, • • • can be 
obtained using Mathematica 9. 



Using 5 0 (x, t) = 1 — £ e~ x as initial approximation. 
We have the solution to problem (1) given as see [] for 
details 



S l (x,t) = h%e- x t-2h% 2 e- 2x t 
S 2 (x,t) = 
h(n + h)^ x t - 2h(n + /*)£ V z '7 - 



(19) 



2_*. h 2 $e- x t 2 



+6h 2 fe- 2x t 2 -9h 2 ^e- 3x t 2 (20) 
S 3 (x, t) = h(n + hf %e~ x t - 
-2h(n 2 -h 2 )% 2 e- 2x t-h 2 (n + h)%e- x t 2 
+I2h 2 (n + h)fe' 2x t 2 - 

-\%h 2 (n + h)^e- 3x t 2 + h ^ e ^ 

6 

28/z 3 £V 2 ¥ C1 , 3e3 -3, 3 176/z 3 £W 

- + 5lh t e t - . 

3 ■ 3 

(21) 

Similarly, S m (x,t) for m = 4, 5, • • • can be 
obtained using MATLAB or Mathematical. 

Then the series solution expression by q-HAM can 
be written in the form 



M 

S(x,t; n; h) = S M (x, t; n;h) = ^ S t (x, t; n; h) 



7=0 



(-] 
(22) 

Equation (22) is an appropriate solution to the 
problem (1) in terms of convergence parameter h and 
n. 

4.3 Numerical Analysis 

4.3.1 Comparison Result of VIM and q-HAM 
In this subsection, we present graphs of solution to 
equation (1) obtained by VIM and q-HAM using 
Mathematica 8 with appropriate h = —0.6 and 
n = 1 . We use only the 2-term approximate solution 
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for both VIM and q-HAM. 

4.3.2 Comparison of q-HAM, VIM and ADM 
We have compared the analytical solution of 
equation (1) obtained by q-HAM, VIM and that of 
ADM in [12] presented in Table 1 and Figs. (3) and 
(4). This is possible when n = 1 and setting h = — 1 .0. 
A good correlation is obtained, so expected 
better convergence if n is larger in the case of 
q-HAM. 

Table 1 



4.3.3 The effect of auxiliary parameter h in the 
solution obtained by q-HAM compared with that of 
VIM 

In this subsection, we give some numerical analysis 
on the effect of auxiliary parameter h in the solution 
obtained by q-HAM compared with the solution 
obtained by VIM. This effect gives it edge over other 
analytical methods of solution to non-linear partial 
differential equations. 




T 


X 


* q-HAM 




$ADM 


0.02 


04000 


0.0145 


0.0144 


0.0146 




0.2000 


5.1 05 1 


0.1091 


0.1052 




0.3000 


0.1377 


0.1877 


0.1878 




04000 


0.2329 


0,2629 


0.2629 




0.5000 


D.3313 


0.3313 


0.3313 




0.6000 


0.3 934 


0.3934 


0.3934 




0.7000 


0.4499 


0.4499 


0.4499 




0.3000 


0.5013 


0.5013 


0.5012 




0.9000 


0,5479 


0.5479 


0,5473 




1.0000 


0,5902 


0.5902 


0.5302 


0.04 


0.1000 


0,0304 


0.0301 


0.0315 




0.2000 


0.1173 


0.1171 


0.1178 




0.3000 


0.1937 


0.1966 


0.1969 




0.4000 


0.2693 


0.2E91 


0.2692 




0.5000 


0.3355 


0.3354 


0.3354 




0.6000 


0.3960 


0.3959 


0.3958 




0.7000 


0.45II 


0.4511 


0.4510 




0.3000 


0,5014 


0.5014 


0.5013 




0.9000 


0,5472 


0.5472 


0.5471 




1.0000 


0.5B90 


0.5883 


0.5888 




Fig. 1 Solution plot obtained by VIM for 0 < t < 0. 1 

and 0 < X < 1 



Fig. 2 Solution plot obtained by q-HAM for 0 < t < 0. 1 

and 0 < X < 1 
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S(HAM) 

S(VIM) 

S(ADM) 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.S 0.3 1 



Fig. 3 The Plot of q-HAM, VIM and ADM for fixed t = 0.02 and 0 < X < 1 




S(HAM) 

S[VIM) 

SfADM} 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Fig. 4 The Plot of q-HAM, VIM and ADM for fixed t = 0.04 and 0 < X < 1 

0.6 h 




Fig. 5 q-HAM solution and VIM solution for fixed t = 0.02, h = -l.5 and 0 < X < 1 
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Fig. 6 q-HAM solution and VIM solution for fixed t = 0.02, h = 1 .5 and 0 < X < 1 




Fig. 7 q-HAM solution and VIM solution for fixed t = 0.02, h = -1 .0 and 0 < X < 1 




Fig. 8 q-HAM solution and VIM solution for fixed t = 0.02, h = 1 .0 and 0 < X < 1 
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Fig. 11 q-HAM solution and VIM solution for fixed t = 0.1, h- -1.5 and 0 < X < 1 
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Fig. 14 



q-HAM solution 



and VIM solution for fixed f = 0.1, h=0.5 and 0<X<1. 
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5. Conclusion 

This paper has successfully developed VIM and 
q-HAM to obtain approximate solution of the nonlinear 
partial differential equation arising from 
fingero-imbibition phenomena in double phase flow 
through porous media. The performance of HAM and 
other method like ADM is greatly improved by q-HAM 
shown in figures presented. We have from our results 

j^J enables the 

q-HAM to match comparatively with that of VIM. The 
efficiency and accuracy of these methods is obvious 
from the graphs, when compared with the result 
obtained in [6] using ADM. 
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Abstract: In the present paper, we introduce the notion of slant submanifolds of an almost hyperbolic contact metric manifolds. We 
have obtained some results on slant submanifolds of an almost hyperbolic contact metric manifolds. We have given a necessary and 
sufficient condition for a slant submanifold of an almost hyperbolic contact metric manifolds. 

Keywords: Almost hyperbolic contact metric manifolds, slant submanifold 

1. Introduction 

Almost contact metric manifold with an almost 
contact metric structure is defined and studied by 
Blairfl]. On the other hand, almost hyperbolic 
,g,r/,^ -structure was given and studied by 
Upadhyay and Dube [7] 

The notion of a slant submanifold of an almost 
Hermitian manifold was given by B.Y.Chen [4]. A. 
Lotta who has defined and studied slant submanifolds 
of an almost contact metric manifold [6]. Later, L. 
Cabrerizo and others investigated slant submanifolds 
of a Sasakian manifolds [3]. The same concept was 
studied under the name "Slant submanifolds of a 
Kenmotsu Manifold" by R.S.Gupta, S.M.K.Haider and 
Mohd.H.Shahid [5].. In 2013, M.Ahmad and K.Ali 
studied "Semi-invariant submanifold of nearly 
hyperbolic cosymplectic manifold" and 
"CR-submanifold of a nearly hyperbolic cosymplectic 
manifold"[2]. 

In this paper we study slant submanifolds of an 
almost hyperbolic contact metric manifolds manifold 
and we find some interesting results. 
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2. Preliminaries 

Let M be a n-dimensional almost hyperbolic 
contact metric manifold with the almost hyperbolic 
contact metric {J>, ?7, g ) -structure, where a tensor <j) 
of type a vector field % and r/ the dual 1-form 

of E, satisfying the followings: 

(/> 2 X = X + ri(X)i; (1) 

g{X,$) = V {X), rj{^-\ (2) 
tf = 0, 770^ = 0 (3) 
g(<?>X,<?>Y) = -g(X,Y)-J 7 (X)j 7 (Y) (4) 

g(0X,Y) = -g(X,0Y) (5) 

for any vector field X and Y in TM . 

Let M be a submanifold immersed in M We 
assume that the vector field E, is tangent to M For 

any X e T (TM) and TV e F (TM L ) , we write 
0X = PX + FX 

0N = tN + fN (6) 

where PX(resp.FX) denotes the tangential (resp. 
normal) component of (f>X , and tN(resp. fN) 
denotes the tangential (resp.normal) component of 
(j)N. If we denote by D the orthogonal distribution to 
E, in TM, we can consider the orthogonal direct 
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decomposition TM = D®{^). 

For each non-zero X tangent to M at X such that 

X isn't proportional to % x , we denote by the 
Wirtinger angle of X that is, the angle between (j>X 
and T x M. Then, M is called slant if is a 

constant, which is independent of the choice of x e M 

a ndXsT x M-{l}[4\ 

Let V be the Riemannian connection on M Then 

V x Y = V x Y + h(X,Y) (7) 

and 



-A N X + V X N 



(8) 



for X,YsT(TM) and NeT^M 1 ]^ and A N 
are the second fundamental form and weingarten 
operator respectively, h and ^4 w are related by 

g(A N X,Y) = g(h(X,Y),N) (9) 
for X,YsT(TM) and N eT[TM L \ 



3. Slant Submanifolds of Almost Hyperbolic 
Contact Metric Manifolds 

In the present section we prove a characterization 
theorem for slant submanifolds of an almost hyperbolic 
contact metric manifolds. 

Theorem 1. Let M be a submanifold of an almost 
hyperbolic contact metric manifold M such that 
^er(rM). Then, M is slant if and only if there 

exist a constant X e [0, l] such that 

P 2 = X(l + r]®<5) (10) 

Furthermore, if 9 is the slant angle of M , then 

X = cos 2 9. 
Proof: : 

=>) : Let M be a slant submanifold of an almost 
hyperbolic contact metric manifold M with slant 
angle 9. % sT(TM). Then, from (6) 

g(PX,</>X) 



= g(PX,PX) 

\\ px \\¥ x \\ 
_ IM 



= const. 



for X G r (TM) . In this case, from (l l) , we have, 
\\PX\\ = \tyX\\. cos 9 (12) 

If we take the square (12) , we get 

g(PX,PX) = g(</,X,</,X).cos 2 9 (13) 
Moreover, 

g(PXJX) 
cos9= ,, V „,, ,/ (14) 

ll PX llll^ X ll 
g(X,^X) 

IMIII^II 

g(X,FPX) g(X,PPX) 
||px||||^x|| ||PX||||^X|| 

g(x,p 2 x) 



cos 9 = 



il^llll^ll 



(11) 



\\ px \¥ x \\ 

all vector field X. Th en from (2),(4),(l2) and 
(14) we find 

cos 0. ||PX|| \$X\ = -g(x, P 2 X) 

cos9.\\</>X\\\\</>X\\.cos9 = -g(x,P 2 x) 

cos 2 9.g(<?>X,<?>X) = -g(x,P 2 x) 

cos 2 9{-g{X,X)-r J {X)r J {X)} = -g(P 2 X,x) 

cos 2 0 {-g (X + 7 (X) £ X)} = -g(P 2 X, X) 

In this case, we have 

P 2 X = cos 2 9{X + T](X)^} 

for any vector field XeT(rM).That is, for all 
XeT(TM) 

P 2 = A(l + rj®%) (15) 
X = cos 2 9 

=>): Let P 2 =X[l + rj<S><^). In this case, from 
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(4) and (10), we find 



COS0 = 



\\ px \\w x \\ 
g(p 2 x,x) 

\\ px \\¥ x \\ 

_ -Zg(X + Tj(X)t,X) 
\\PX\\\\^X\\ 
g(X,X) + rj(X)tj(X) 

\\ px \\¥ x \\ 

_ o g(0X,0X) 

\\ px W x \\ 
¥ x \\ 



= -A 
■■A 
= A 



\PX\ 

for any vector field X e r (TM) . On the other hand, 

ILpxII 

COS# = 



UX\\ 



1 



= A- 

COS& 

then we have 

COS 2 0 = A (16) 
Since A is constant, 6 is constant too. In this case, 
M is slant. 

Corollary 2. Let M be a slant submanifold of an 
almost hyperbolic contact metric manifold M , with 
slant angle 6 and %&T(TM). Then, 
g(PX,PY) = -cos 2 6{g(X,Y) + r 1 (X)r 1 (Y)} 

(17) 

g(FX,FY) = -sin 2 e{g{X,Y) + ri{X)ri{Y)} 

(18) 

for any vector field X and Y in TM. 

Proof: We can find with basic calculations using 
Theorem 3.1. 

Theorem 3. Let M be a submanifold of an almost 
hyperbolic contact metric manifolds M such that 
^eT(TM). in this case, 

VQ = 0 (19) 

Proof: From (lO), we get 



QY = cos 2 0{Y + ?](Y)^} (20) 

for any vector field Y in TM. If we take the 
covariant derivative of (17) along X , we get 

V x (er) = cos 2 #{v x (r + /7(rK)} (21) 

= cos 2 0 { v x F + v x (ri ( Y) % + rj ( Y) yj)} 

= cos 2 0{v x Y + ?](v x Y)t} 
for X, Y e Y (TM ) . On the other hand, 

[v x Q)y = vAQ y )-QvxY 

= cos 2 #{v;/ + /7(v^K}- 

- cos 2 £{^ + 77(^)^=0 
for any vector fields X,Y in r(77kf). 

Theorem 4. Let M be a submanifold of an almost 
hyperbolic contact metric manifolds M such that 
•Then, the endomorphism Q \d has only 
one eigenvalue A, at each point of M Moreover 
A = cos 2 6. 
Proof: From (10), 

QX = P 2 X = cos 2 0{X + T](X)^} (22) 

= (cos 2 e)x 
= AX 

for all vector field X in TM. From (22), we see that 

the endomorphism Q\ D has only one eigenvalue. 
Theorem 5. Let M be a submanifold of an almost 
hyperbolic contact metric manifolds M such that 
£ e T (TM) .M i s slant if and only if 

a-) The endomorphism 2 Id has only one 
eigenvalue at each point of M . 

b-) (vx<2K = o 

for any vector fields X and F . 

Proof: =>) : Suppose that M is slant. Then, 

statements (a) and (b) follow directly from 

Theorem 3.2 and Theorem 3.3. 

<=) : Let Ay (x) be the eigenvalue of Q \ D at each 
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point x of M and Y e D be a unit eigenvector 
associated with \, i.e., 

QY = \Y (23) 

If we take the covariant derivative of (23) along 
X , we find 

V X {QY) = V X [\)Y + \V X Y (24) 

for any vector field X ,Y in TM. From (b), we get 

(vxQ)y=vAQ y )-Qv* y 

= 0 

Vx{QY) = Qv x Y (25) 
all X,YsT(TM). From (24) an d (25), we find 

QVxY = VxMY + \v x Y (26) 

all X,YeT(TM). Since both VV 7 an d Qv X Y 
are perpendicular to Y, from (26), we get 

X(4) = 0 (27) 

for all X e r(JM). Then \ is constant on M . 

Now we prove that M is slant. In view of (13), it 
is enough to show that there exists a constant ju such 

that 

0 = A#(/ + i7®#) 

Let X be in TM. Then 

X = X-77(X)£e£> 

Hence 



Since 

Q\ D =\i 

we have 
and so 

QX = \X = X l (X + i 1 (X)$) 

Theorem 6. By taking ju = \ , we get that M is 
slant. Moreover, if M is slant, X - \ - JU - COS 2 6, 
where 9 denotes the slant angle of M . 
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